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ABSTRACT 


Seismic  studies  of  the  oceanic  crust,  both  experimental  and  theoretical,  often 
assume  a  flat  seafloor  and  laterally  homogeneous  crust.  This  is  done  regardless  of  the 
appearance  in  seismic  data  of  obvious  effects  due  to  scattering  from  lateral  heterogeneities 
both  on  and  in  the  seafloor.  Detailed  fine  scale  surveys  of  mid-ocean  ridges,  where  the 
upper  oceanic  crust  is  exposed,  have  revealed  the  presence  of  lateral  heterogeneities  in  the 
form  of  complicated  topography,  extrusive  volcanic  structure,  and  abundant  fracturing  and 
faulting.  These  heterogeneities  have  a  significant  affect  on  the  propagation  of 
seismo/acoustic  energy  through  the  crust,  especially  in  the  immediate  vicinity  of  the 
seafloor.  This  thesis  deals  with  the  problem  of  scattering  of  seismo/acoustic  energy  from  a 
number  of  forms  of  lateral  heterogeneity  in  the  upper  oceanic  crust. 

A  common  theme  throughout  this  work  is  that  the  size  of  the  heterogeneity  on  or  in 
the  seafloor  is  of  the  same  order  of  magnitude  as  the  seismo/acoustic  wavelength.  This  is 
the  realm  of  scattering  theory  where  the  wave-like  characteristics  of  seismic  energy  have  a 
particularly  large  influence  on  the  outcome  of  interaction  with  structure  in  the  media.  The 
work  presented  here  involves  the  application  of  the  finite  difference  modeling  technique  to 
problems  concerning  laterally  heterogeneous  elastic  media.  This  method  is  a  full  wave 
solution  to  the  elastic  wave  equation  and  as  such  includes  all  wave  interactions  with  the 
media.  The  finite  difference  formulation  is  used  to  study  four  distinct  phenomena; 
scattering  from  discrete  deterministic  seafloor  features;  wave  propagation  through 
continuous  randomly  heterogeneous  upper  oceanic  crust;  scattering  from  more  complicated 
topographic  profiles  and  the  limitations  of  the  method  for  the  rough  seafloor  problem;  and 
the  problem  of  plane  acoustic  wave  scattering  from  an  infinite  elastic  cylinder. 

The  principal  finding  of  this  work  is  that  lateral  heterogeneities  in  the  upper  oceanic 
crust  can  have  a  dramatic  affect  on  seismo/acoustic  wave  propagation.  Scattering  from 
rough  seafloors  and/or  volume  heterogeneities  is  often  quite  similar  and  causes  the 
occurrence  of  signal  generated  'noise'  (coda),  decorrelation  of  primary  arrivals,  and 
anomalies  in  arrival  travel  time  and  amplitude.  Topographic  and  volume  scatterers  acting  as 
secondary  sources  of  seismic  energy  can  cause  a  resonant  coupling  of  body  wave  energy 
into  interface  (Stoneley)  waves  at  the  seafloor.  This  is  possibly  one  mechanism  by  which 
natural  seismic  and  storm  generated  acoustic  energy  can  be  coupled  into  seafloor  noise. 
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The  applicability  of  the  use  of  the  finite  difference  method  for  non-planar  water- 
solid  interfaces  is  also  discussed.  Models  were  calculated  which  approximate  sinusoidal 
seafloors  and  plane  acoustic  wave  scattering  from  an  infinite  elastic  cylinder.  The 
discretization  of  a  rectangular  difference  grid  must  be  extremely  fine  to  accurately 
accommodate  a  smoothly  varying  water-solid  interface  which  does  not  align  with  the  grid. 
Regardless  of  the  discretization  concerns,  the  rough  seafloor  models  presented  here 
demonstrate  the  arrivals  expected  from  larger  scale  sinusoidal  topography  as  well  as  the 
importance  of  considering  quite  small  (<1/15  wavelength)  topographic  features  in  the 
scattering  problem.  Also,  steep  topography  will  allow  seismo/acoustic  energy  to  enter  the 
seafloor  at  very  large  ranges  because  the  angle  of  incidence  can  repeatedly  fall  below  the 
critical  angle  for  transmitted  energy,  especially  for  convened  shear  energy.  Ray  theoretical 
shadow  zones  do  not  occur  in  these  models  (or  in  the  real  world)  because  of  Franz-type 
waves  diffracting  into  areas  where  the  grazing  angle  is  less  than  zero. 
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The  ultimate  objective  of  much  of  the  work  in  the  fields  of  marine  geology  and 
geophysics  is  to  understand  the  structure  and  processes  of  evolution  of  the  oceanic  crust. 
Marine  seismology  is  the  science  of  the  use  of  sound  wave  energy  (both  acousnc  and 
seismic)  to  remotely  examine  the  structure  of  the  ocean  subbottom.  While  the  exact  nature 
of  the  relationship  between  seismic  velocity,  impedance,  and  lithology  is  still  under 
investigation  (Karson  and  Fox,  1986;  Lewis,  1983;  Spudich  and  Orcutt.  1980:  Spudich 
and  Orcutt,  1980) ,  there  is  no  doubt  that  the  results  of  marine  seismic  experiments  have 
told  us  much  about  the  gross  structure  of  the  oceanic  crust  Inherent  in  the  use  of 
seismology  as  a  tool  is  the  assumption  of  a  basic  understanding  of  the  physics  of  wave 
propagation  through  the  medium  in  question.  The  physics  of  seismic  wave  propagation 
through  elastic  media  have  been  used  in  various  states  of  approximation  to  model  seismic 
data  collected  in  the  field.  This  thesis  involves  the  use  of  the  finite  difference  modeling 
technique  to  study  wave  propagation  through  laterally  heterogeneous  upper  oceanic  crust. 

The  classical  picture  of  the  oceanic  crust  is  one  of  a  discrete  number  of  flat-lying 
homogeneous  layers  with  fairly  predictable  velocities  which  correlate  with  age  (Houtz  and 
Ewing,  1976).  This  layered  structure  is  generally  composed,  in  descending  order,  of 
sediments,  extrusive  basalts,  intrusive  gabbroic  dikes,  and  a  peridotitic  mantle  at  a  depth  of 
5-7  kilometers  below  the  sediment-basement  interface.  Simple  seismic  models  using 
velocities  assigned  to  these  layers  were  sufficient  to  reproduce  primary  wave  travel  times  of 
early  seismic  experiments.  More  accurate  modeling  techniques  and  better  controlled  marine 
seismic  experiments  have  since  shown  that  the  classical  oceanic  structure,  while  close  to 


correct  'on  average',  does  not  explain  the  lack  of  inter-crustal  reflections  and  the  amplitude 
of  refracted  arrivals  seen  in  the  field  data  (Spudich  and  Orcutt,  1980).  Rather,  the  'layers’ 
of  oceanic  crust  are  now  delineated  as  areas  of  crust  with  fairly  constant  vertical  velocity 
gradients  (Collins,  1988;  Spudich  and  Orcutt,  1980).  More  sophisticated  modelling 
techniques  can  now  reproduce  travel  times,  amplitudes,  and  phase  information  seen  in 
seismic  data  (Helmberger  and  Morris,  1970;  Fuchs  and  Muller,  1971). 

Seismo/acoustic  modelling  and  state-of-the-art  seismic  experimental  methods  have 
traditionally  driven  each  other  to  higher  levels  of  accuracy.  Although  the  vertical  picture 
has  changed,  lateral  homogeneity  in  the  upper  oceanic  crust  is  still  assumed  today  in  most 
cases.  This  is  due  in  part  to  the  ability  of  1 -dimensional  velocity  structures  to  reproduce  the 
gross  features  of  seismograms,  but  is  also  due  to  the  inability  of  most  '2-dimensional' 
modelling  schemes  to  handle  truly  2-dimensional,  or  laterally  heterogeneous,  earth  models. 
There  are  secondary  features  in  virtually  all  seismic  data  which  are  due  to  scattering  from 
some  type  of  lateral  velocity  feature. 

Another  commonly  held  belief  is  that  seismic  methods  cannot  resolve  structure 
which  has  lateral  extent  on  the  order  of  the  seismic  wavelength  in  size.  This  is  true  in  the 
context  of  trying  to  resolve  laterally  homogeneous  structure  with  vertical  extent  less  that  a 
wavelength  in  size.  However,  heterogeneities  of  this  size  have  a  dramatic  effect  on  seismic 
wave  propagation.  In  fact,  scattering  of  primary  energy  reaches  a  maximum  when  ka( k  = 
wavenumber  =  2k  /  wavelength,  a=heterogeneity  length  scale)  is  approximately  equal  to 
one  (Aid  and  Richards,  1980).  Scattering  of  energy  from  heterogeneities  manifests  itself 
mainly  in  secondary  features  of  the  seismograms  such  as  signal  generated  'noise'  (coda), 
secondary  body  and  interface  waves,  and  signal  decorrelation.  If  the  heterogeneities  are 
large  enough,  deterministic,  predictable  effects  such  as  noticeable  travel  time  and  amplitude 
anomalies  will  also  appear  in  the  seismograms. 
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This  dissertation  is  a  study  of  the  effects  of  a  laterally  heterogeneous  upper  oceanic 


crust  on  the  scattering  of  seismo/acousnc  energy  both  at  and  below  the  seafloor.  The 
method  of  finite  differences  is  used  to  model  propagation  effects  in  media  with 
heterogeneity  sizes  on  the  same  order  of  magnitude  as  the  seismic  wavelength.  Chapters  in 
this  work  cover  a  number  of  types’  of  lateral  heterogeneity.  Scattering  ffom  deterministic, 
isolated,  seafloor  features  is  investigated  in  Chapter  1 .  Chapter  one  appears  here  with 
permission  from  Journal  of  the  Acoustical  Society  of  America  (Dougherty  and  Stephen, 
1987).  Crustal  volume  heterogeneities  can  also  have  a  strong  influence  on  propagation. 
Chapter  2  deals  with  these  types  of  media  and  also  appears  here  with  permission  ffom  the 
Journal  of  Pure  and  Applied  Geophysics  (Dougherty  and  Stephen,  1988).  Further 
investigation  into  the  effects  of  more  complicated  topography  is  made  in  chapter  3. 

Chapter  4  is  a  comparison  between  the  analytical  and  finite  difference  solutions  to  the 
problem  of  plane  acoustic  wave  scattering  ffom  an  infinite  elastic  cylinder. 

Chapter  one  addresses  the  problem  of  deterministic  scattering  from  distinct  seafloor 
features.  Seismic  data  often  contain  sections  of  relatively  large  anomalous  arrivals 
superimposed  on  expected  arrivals  for  small  groups  of  traces.  One  line  of  airgun  data  ffom 
an  ocean  bottom  hydrophone  (OBH)  in  the  Rivera  Ocean  Seismic  Experiment  (ROSE) 
contains  a  very  coherent  diffracted  arrival  just  after  the  primary  refracted  wave.  This  arrival 
is  referred  to  as  a  "refraction  branch  diffraction"  to  distinguish  it  ffom  diffraction 
hyperbola  of  near  normal  incidence  reflections.  The  finite  difference  method  (Bhasavanija, 
1983;  Clayton  and  Engquist,  1977;  Kelly  et  al.,  1976;  Nicoletis,  1981;  Stephen,  1984)  was 
used  to  demonstrate  that  there  are  no  simple  reflection  paths  which  can  account  for  the 
arrival.  Rather,  it  is  due  to  secondary  scattering  from  any  of  a  number  of  seafloor  features 
with  height  about  equal  to  the  acoustic  wavelength.  Energy  is  scattered  ffom  the  seafloor 
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features  both  into  the  water  column  and  back  into  the  crust  where  it  propagates  to  the 
receiver. 

The  interaction  of  seismo/acousdc  energy  with  discrete  scatterers  is  fairly  simple 
conceptually  when  compared  to  the  problem  of  propagation  through  a  continuously  varying 
medium.  Random  distributions  of  heterogeneity  in  the  upper  oceanic  crust  can  account  for 
much  of  the  signal  generated  'noise'  often  seen  in  marine  seismic  data.  Also,  decorrelation 
in  the  primary  arrival  waveforms  is  caused  by  interaction  with  strongly  heterogeneous 
media.  These  effects  are  studied  in  chapter  2  by  calculating  the  wave  propagation  through  a 
number  of  marine  models  with  random  velocity  perturbations  containing  a  range  of 
heterogeneity  correlation  lengths. 

The  theoretical  study  of  wave  propagation  through  media  with  continuous  volume 
heterogeneities  has  been  well  developed.  However,  the  solution  to  the  problem  of  the 
influence  of  strong  velocity  variations  with  ka  near  one  remains  to  be  a  difficult  one. 
Historically,  most  of  the  effort  concerning  scattering  from  volume  heterogeneities  has  been 
directed  at  the  amplitude  attenuation  and  coda  wave  excitation  of  seismic  energy  from 
earthquakes  or  teleseismic  explosive  sources  (Aki,  1973;  Aki,  1982;  Frankel  and  Clayton, 
1986;  McLaughlin,  Anderson  et  al.,  1985;  McLaughlin,  Johnson  et  al.,  1983;  Menke  and 
Chen,  1984;  Menke  and  Dubendorff,  1985;  Sato,  1982;  Wu,  1982;  Wu,  1983;  Wu  and 
Aki,  1988).  Large  scale  heterogeneities  have  the  greatest  influence  on  primary  arrival 
amplitude  attenuation  while  smaller  scale  heterogeneities  seem  to  control  coda  wave 
frequency  content  and  decay  (Aki,  1973;  Aki,  1982).  Most  analytical  work  done  to  explain 
these  results  depends  on  the  validity  of  the  Bom  approximation  which  only  considers  the 
first  interaction  of  energy  with  a  scatterer  (no  multiple  scattering)  and  that  the  scattering  is 
weak  with  respect  to  the  primary  wave  amplitude.  When  the  experiment  is  scaled  to  that  of 
typical  marine  refraction  work,  the  Bom  approximation  is  not  necessarily  valid  (Hudson, 
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1982;  Hudson  and  Heritage,  1981).  The  case  of  scattering  of  vector  waves  from  a  point 
source  travelling  through  a  medium  with  velocity  gradients,  topography,  and  random 
heterogeneities  has  little  chance  to  be  solved  analytically.  Therefore,  numerical  methods  are 
probably  the  imst  effective  resource  for  studying  this  problem..  The  finite  difference 
method  has  already  been  applied  to  the  problem  for  the  case  of  laterally  heterogeneous 
halfspaces  with  plane  wave  sources  (Frankel  and  Clayton,  1986;  McLaughlin  et  al.,  1985; 
McLaughlin  et  al.,  1983).  However,  these  studies  have  not  dealt  with  a  pulse  source 
travelling  through  the  strong  heterogeneities  or  the  water-solid  interface  of  the  marine  case. 
The  marine  case  is  particularly  interesting  because  of  the  coupling  of  acoustic  and  seismic 
energy  at  the  seafloor.  Stephen  has  used  this  method  successfully  to  study  the  scattering 
effects  of  simple  later  J  heterogeneity  in  the  oceanic  crust  (Stephen,  1984a,b;  1985; 
1988a,b,c;  Stephen  and  Bolmer,  1985).  The  finite  difference  formulation  used  for  the 
work  in  chapters  2  and  3  is  a  displacement-stress  formulation  and  is  well  suited  for  use 
with  the  marine  problem  (Stephen,  1988;  Virieux,  1984;  Virieux,  1986). 

In  general,  random  scattered  'noise'  in  the  seismograms  and  decorrelation  of  the 
primary  compressional  diving  wave  reaches  a  maximum  for  ka  near  one.  As  ka  increases 
much  above  one,  deterministic  travel  time  and  amplitude  effects  become  more  prevalent. 

As  ka  decreases  much  below  one,  the  medium  becomes  an  equivalent  homogeneous 
medium  (Aki  and  Richards,  1980)  and  coherence  of  primary  arrivals  returns.  Another 
interesting  result  of  chapter  2  is  that  the  heterogeneities  near  the  water-solid  interface  can  act 
as  secondary  sources  of  interface  waves  when  energy  is  scattered  from  them.  It  has  been 
proposed  that  seafloor  noise  propagates  as  some  type  of  interface  wave  along  the  seafloor. 
However,  the  source  for  these  interface  waves  was  unknown.  The  work  in  chapter  2 
implies  that  this  'noise'  can  be  generated  by  secondary  scattering  of  acoustic  waves  or 
microseisms  from  heterogeneities  within  the  crust  but  near  the  water-solid  interface. 
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Chapter  3  involves  more  investigation  of  scattering  form  rough  seafloors.  The 
problem  of  seismic  energy  scattering  from  heterogeneities  within  the  crust  is  complicated 
by  the  presence  of  topography  on  the  seafloor  or  sediment/basement  interface.  While 
heterogeneities  within  the  crust  can  contain  velocity  contrasts  on  the  order  of  1-20  per  cent, 
the  water-solid  interface  can  represent  a  velocity  contrast  of  over  100  per  cent.  Because  of 
this  sharp  impedance  contrast,  non-planar  seafloors  can  cause  large  scattering  effects  when 
energy  propagates  into  the  oceanic  crust.  With  increased  interest  in  the  seismic  study  of 
mid-ocean  ridges  where  the  impedance  contrast  is  particulaily  high,  it  it  important  to 
understand  the  relationship  between  topography  and  seismo/acoustic  propagation  through 
the  crust. 

As  with  the  problem  of  wave  propagation  through  volume  heterogeneities,  the 
rough  surface  scattering  problem  has  also  been  very  well  studied.  However,  the 
approximations  used  for  most  of  the  known  solutions  make  them  inappropriate  for  the 
study  of  elastic  scattering  at  realistic  seafloors  (Ogilvy,  1987).  These  solutions  include 
Rayleigh's  method  of  plane  wave  summation  (Rayleigh,  1878),  the  method  of  small 
perturbations  (Gilbert  and  Knopoff,  1960;  very  small  topography  with  respect  to  acoustic 
wavelength),  the  Kirchhoff  method  (Eckart,  1953;  assumes  scattering  is  from  tangent 
planes  along  the  topography)  and  methods  which  sum  contributions  from  uistributions  of 
point  sources,  facets,  or  protuberances  along  the  bottom  (see  Ogilvy,  1987). 

Methods  which  require  fewer  limiting  assumptions  and  are  more  flexible  have  more 
recently  been  used  to  study  this  problem.  Bouchon  (1985),  Campillo  and  Bouchon 
(1985),  and  Paul  and  Campillo  (1988)  have  successfully  used  a  boundary  integral 
formulation  coupled  with  the  discrete  wavenumber  method  to  study  diffractions  from 
corrugated  boundari  ;s.  The  computationally  intensive  finite  difference/element,  and 
pseudospectral  method  (Kosloff,  et  al.,  1984;  Fomberg,  1986,  1988)  have  also  become 
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attractive  because  of  the  full  solution  to  the  problem  of  surface  and  subsurface  scattering 
with  very  few  initial  limiting  assumptions.  Stephen  (1984)  and  the  work  from  chapter  1  of 
this  work  use  the  method  of  finite  differences  to  explain  scattering  effects  of  discrete 
seafloor  features  on  both  sides  of  the  critical  range.  Hill  and  Levander,  (1984)  and 
Levander  and  Hill  (1985)  also  used  the  finite  difference  method  to  model  scattering  from 
rough  low  velocity  layers  within  the  crust  and  rough  elastic  surface  layers. 

In  chapter  3  the  finite  difference  method  is  used  to  model  propagation  through 
sinusoidal  seafloors.  While  the  initial  intent  of  this  work  was  to  help  to  distinguish 
between  volume  and  surface  scattering,  this  chapter  actually  deals  mainly  with  the  question 
of  the  appropriateness  of  the  use  of  the  method  for  this  type  of  problem.  The  principle 
problem  with  these  models  is  in  the  approximation  of  a  curved  seafloor  surface  on  a 
rectangular  grid.  A  rectangular  stair-stepped  structure  must  be  used  to  closely  approximate 
a  continuously  sloping  seafloor.  In  this  way,  a  kind  of  secondary  microroughness  is 
imposed  onto  the  larger  scale  sinusoidal  topography.  Even  though  the  stair  steps  are  quite 
small  in  relation  to  the  seismic  wavelength,  significant  diffractions  still  occur  from  the 
microroughness  on  the  seafloor.  Scattering  phenomena  due  to  this  problem  are  discussed 
in  chapter  3  as  well  as  some  'real'  effects  which  can  be  expected  from  sinusoidal 
topography. 

The  problem  of  plane  acoustic  wave  scattering  from  an  infinite  elastic  cylinder  is 
solved  in  chapter  4  with  the  finite  difference  formulation  used  in  the  work  of  chapters  2  and 
3.  The  intent  of  this  work  was  not  to  shed  any  new  light  on  the  solution  to  the  problem, 
but  rather,  to  use  this  well  studied  problem  with  exact  analytical  solutions  as  a  rigorous  test 
of  the  method  and  formulation.  Work  on  the  cylinder  problem  pointed  out  a  slight 
asymmetry  in  the  coding  of  the  definition  of  elastic  parameters  on  the  finite  difference  grid 
for  the  displacement-stress  formulation.  The  coding  mistake  was  corrected  and  symmetry 

15 


of  the  solution  was  established.  All  of  the  models  in  chapter  2  and  the  preliminary  models 
of  chapter  3  were  calculated  before  this  slight  asymmetry  in  the  grid  was  discovered. 
Fortunately,  the  effect  is  a  very  small  one  and  test  reruns  of  some  of  these  models  with  the 
corrected  grid  definition  showed  no  measurable  difference  in  mode!  results  discussed  in  the 
text.  The  TEST  models  presented  in  chapter  3  and  all  of  the  models  presented  in  chapter  4 
were  run  with  the  corrected  grid  definition. 
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ABSTRACT 

A  strong  "refraction  branch  diffraction",  presumably  due  to  scattering  from  a  lateral 
heterogeneity  on  or  below  the  sea  floor,  has  been  observed  on  ocean  bottom  hydrophone 
data  from  the  Rivera  Ocean  Seismic  Experiment  (ROSE).  This  arrival  is  unusual  because 
of  its  coherence  and  relatively  large  amplitude.  Finite  difference  modeling  of  a  number  of 
possible  seafloor  diffractors  and  associated  lateral  velocity  variations  are  presented  which 
demonstrate  the  occurrence  and  characteristics  of  "refraction  branch  diffractions".  In 
general,  the  half-width  of  the  diffractor  must  be  approximately  the  same  as  the  seismic 
wavelength  in  order  to  produce  a  strong  diffraction.  Velocity  gradients  present  in  the 
models,  as  well  as  P-S  conversion,  complicate  the  wavelength-halfwidth  relationship. 
Three  different  models,  a  hill,  a  valley,  and  a  subsurface,  high-velocity  block,  all 
produced  diffractions  of  sufficient  amplitude  to  explain  the  data.  There  is  a  hill  along  the 
line  with  approximately  the  same  dimensions  as  the  model  hill  and  it  is  the  proposed 
source  of  the  diffracted  energy  in  the  data.  The  large  models  used  also  clearly  demonstrate 
the  existence  of  phases  which  are  theoretically  possible  but  rarely  seen  in  marine  seismic 
(geoacoustic)  data  such  as  the  pseudo-Rayleigh  wave  and  the  P  and  S  interference  head 
waves. 
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INTRODUCTION 


Marine  seismic  (geoacoustic)  data  often  contain  relatively  large  anomalous  arrivals 
superimposed  on  expected  arrivals  for  individual  or  small  groups  of  traces.  Factors  such 
as  topographic  focusing  or  crustal  heterogeneities  could  be  the  cause  of  such  anomalies. 
Ocean  bottom  hydrophone  (OBH)  data  from  the  Rivera  Ocean  Seismic  Experiment 
(ROSE)  contain  much  of  this  type  of  noise.  There  is  one  OBH  line,  however,  which 
contains  coherent  arrivals  between  the  refracted  P-wave  and  the  direct  water  wave  for 
almost  the  entire  ten  kilometer  section  (figure  1.1).  We  call  this  the  "refraction  branch 
diffraction”  to  distinguish  it  from  diffraction  hyperbola  of  near  normal  incidence 
reflections.  There  are  no  simple  reflection  paths  which  can  account  for  the  travel  times  of 
the  refraction  branch  diffraction.  The  arrival  time  behavior  indicates  that  some  type  of 
diffraction  or  back-scattering  is  occurring  at  a  range  of  6.6  kilometers.  Rough  seafloor  in 
the  area,  seen  in  figure  1.2,  contains  many  hills  and  valleys,  one  of  which  is  the  most 
likely  cause  of  this  scattering.  In  this  paper  we  model  and  discuss  the  phenomenon  of  the 
refraction  branch  diffractions. 

The  ability  to  detect  heterogeneities  in  or  on  the  oceanic  crust  is  a  function  of  the 
frequency  of  the  seismic  source  (or  the  wavenumber  or  wavelength  of  the  energy  as  it 
travels  through  the  medium),  the  length  scale  of  the  heterogeneity  and  the  distance  travelled 
by  the  energy  through  the  heterogeneous  material.  Aki  and  Richards  1  classify  scattering 
phenomena  on  the  basis  of  the  dimensionless  parameters  ka  (2tc  times  the  ratio  of 
heterogeneity  length  to  wavelength)  and  kL  (2k  times  the  number  of  wavelengths  travelled 
through  the  heterogeneous  region).  Small  scale  heterogeneities  (small  ka)  can  be 
approximated  by  an  equivalent  homogeneous  bodyl  and  the  scattering  can  be  considered  as 
an  attenuation  effect.  Large  scale  heterogeneities  (large  ka)  can  be  considered  as  separate 
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homogeneous  bodies  and  ray  theory  can  be  used  to  trace  the  energy  propagation. 

However,  when  the  size  of  the  heterogeneity  is  approximately  equal  to  the  seismic 
wavelength,  energy  loss  due  to  scattering  reaches  a  maximum.  For  the  problem 
considered  here  (figure  1.1)  it  appears  that  there  is  a  single  scatterer  ( kL  small,  lea  near  1 ) 
and  modeling  in  a  deterministic  fashion  (eg.  using  the  method  of  Finite  differences)  is 
appropriate. 

For  the  purpose  of  travel  time  analysis,  a  scattering  body  can  be  considered  as  a 
secondary  source  and  ray  theories  can  be  used  for  modeling.  These  methods  generally 
assume  either  high  frequencies  or  laterally  homogenous  media  or  both.  Amplitude 
analysis  and  energy  partitioning  in  media  which  varies  laterally  on  the  scale  of  seismic 
wavelengths  are  not  handled  well  with  these  methods.  At  realistic  frequencies  and  with 
laterally  heterogeneous  media,  the  wave  characteristics  of  seismic  propagation  become 
important.  Modeling  schemes  which  preserve  wave  phenomena  such  as  scattering,  energy 
partitioning,  and  interference  must  be  used  under  these  circumstances.  Synthetic  modeling 
using  the  wavenumber  integral  or  finite  difference/element  methods  preserve  wave 
characteristics  but  only  the  finite  difference/element  methods  allow  for  lateral  as  well  as 
vertical  variations  in  structure.  The  finite  difference  method  is  used  in  this  study  to  model 
the  rough  sea  floor  found  in  the  project  ROSE  area. 
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Figure  1.1.  The  refraction  profile  for  OBH  2W  on  airgun  line  1  from  the  ROSE 
experiment  (see  figure  1.3  for  location).  The  bathymetry  section  below  the  seismogram 
was  used  for  the  topographic  travel  time  correction.  The  "refraction  branch  diffraction" 
shows  as  a  "smile"  after  the  primary  refracted  arrival  (modified  figure  from  Purdy4). 
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Figure  1.2.  Echo  sounder  (3.5  kHz'!  record  along  airgun  line  1  near  OBH  2W.  Note  the 
hill  offset  from  the  line  behind  the  large  valley  centered  around  7  km  range.  This  hill  is  the 
proposed  source  of  the  scattered  energy  seen  in  figure  1.1. 
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DATA  AREA 


The  sea  based  portion  of  the  Rivera  Ocean  Seismic  Experiment  took  place  in  early 
1979  and  was  carried  out  to  study  the  evolution  of  young  ocean  lithosphere  at  the  East 
Pacific  Rise  (EPR)  spreading  center  (figure  1.3).  Earthquake  and  explosive  events  from 
the  EPR,  the  Clipperton  and  Orozco  fracture  zones  ,  and  along  the  coast  of  Mexico,  were 
recorded  for  a  two  month  period  using  ocean  bottom  seismometers  (OBS)  and  ocean 
bottom  hydrophones  (OBH)“.  Crustal  ages  of  0-4  million  years  are  covered  by  the  ROSE 
study  area. 

This  study  is  an  attempt  to  model  the  diffraction  seen  in  the  data  from  OBH  2W  on 
airgun  line  1.  Four  7.604  liter  (464  cu.in.)  Bolt  1500C  air  guns  were  used  to  shoot  the 
line,  with  five  OBH's  for  data  recording  (figure  1.3).  The  crustal  age  below  OBH  2W  is 
approximately  0.5  million  yearsA  Data  from  this  line  were  used  primarily  by  Purdy4  to 
determine  the  variability  in  the  seismic  structure  of  oceanic  crustal  layer  two.  Ewing  and 
Purdy^  also  used  part  of  this  data  to  constrain  the  velocity  structure  of  the  upper  500-800 
meters  of  the  oceanic  crust.  Velocity-depth  profiles  used  to  model  the  crust  around  OBH 
2W  are  given  in  figure  1.4.  Separate  upper  crustal  velocity  structures  of  Purdy4  and 
Ewing  and  Purdy^  are  shown  combined  with  a  lower  crustal  P-velocity  of  6.0  km/sec.  S- 
wave  velocities  shown  in  figure  1.4  were  inferred  by  considering  the  upper  crust  as  a 
Poisson  solid  (Poisson's  ratio  (s)  =  0.25).  The  upper  velocity  structure  of  Purdy4  was 
used  for  most  of  the  modeling  in  this  study. 

Figure  1.2  shows  a  portion  of  a  bathymetric  profile  taken  along  airgun  line  1  in  the 
vicinity  of  OBH  2W.  In  general,  the  terrane  around  OBH  2W  is  rough.  The  striking 
feature  of  the  bathymetry  around  OBH  2W  is  the  presence  of  a  relatively  large  valley  with 
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Figure  1.3.  Location  of  OBH  2  along  airgun  line  1  of  the  ROSE  experiment.  Refraction 


data  seen  in  figure  1  (survey  OBH  2W)  includes  shots  ranging  from  0  to  9  km  west  of 
OBH  2.  Inset  shows  general  location  of  the  ROSE  study  area  (modified  figure  from  Ewing 
and  Meyer^). 


31 


VELOCITY  (km/sJ 


Figure  1.4.  Velocity-depth  profiles  used  for  finite  difference  models.  Profile  ’A’  was 
used  for  FLAT,  HILL,  VALLEY1,  VALLEY2,  and  BLOCK  models,  and  profile  B'  was 
used  for  VALLEY3. 
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a  hill  offset  behind  it.  Large  valleys  such  as  this  one  are  generally  not  seen  along  the  rest  of 
the  line.  While  there  are  many  hill  structures  along  the  line,  little  can  be  said  about  their 
three  dimensional  extent,  and  no  initial  conclusions  can  be  drawn  from  the  bathymetry 
alone.  It  is  possible,  for  example,  that  the  dark  hills  in  figure  1.2  are,  in  fact,  lineated 
structures  and  the  hill  offset  behind  the  valley  is  circular.  There  is  some  3-D  seabeam  data 
along  the  crest  of  the  EPR,  but  the  data  stops  jus:  to  the  west  of  the  area  of  OBH  2W 
(Macdonald,  1985,  personal  communication). 

TRAVEL  TIME  ANALYSIS 

The  seismogram  from  OBH  2W  (figure  1. 1)  is  corrected  for  local  ocean  bottom 
topography  using  the  method  of  Purdy 6.  The  diffracted  arrival  appears  symmetric  about  a 
range  of  approximately  6.6  kilometers.  The  minimum  difference  in  travel  time  (DT) 
between  the  refracted  and  diffracted  arrivals  is  about  0.13  seconds  (at  6.6  km)  and  the 
maximum  DT  for  the  ranges  available  is  around  0.53  seconds  (at  4.6  and  8.6  km).  These 
DTs  are  subjective  and  depend  on  the  correct  picks  for  the  incidence  of  the  diffraction.  It 
is  clear,  however,  that  the  shape  of  the  diffraction  is  fairly  well  determined  and  that  only 
its  position  in  time  will  change  with  different  picks.  The  general  shape  of  the  diffracted 
arrival,  as  well  as  the  DTs  were  modeled  first  with  simple  ray  theory.  Also,  the  effects  of 
moving  the  diffracting  body  along  and  out  of  the  sagittal  plane  of  the  model  were  studied. 
These  analyses  follow  for  both  sea  floor  and  subsurface  diffractors. 

For  the  purpose  of  travel  time  analysis,  the  sea  floor  hill  and  the  buried  body  were 
treated  as  point  diffractors  which  act  as  secondary  sources  of  seismic  energy.  That  is, 
when  plane  wave  energy  reaches  the  point  diffractor,  it  is  excited  and  radiates  energy  in  all 
directions.  This  is  opposed  to  refraction  through  the  body  which  would  send  out  energy 
in  only  one  direction,  or  for  a  sharp,  fiat  interface,  in  a  finite  number  of  directions,  from 
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any  one  incoming  plane  wave  front.  Stephen^  used  the  finite  difference  method  to  study 
refractions  through  a  sea  floor  hill.  Although  the  geometry  of  Stephen's  model  is  different 
than  that  used  for  this  study,  it  is  obvious  that  there  would  be  no  back-scattering  of  energy 
from  the  hill  by  considering  only  refraction  phenomena.  Rather,  diffraction  from  a  sea 
floor  hill  is  necessary  to  explain  the  back-scattered  energy.  Size  and  shape  of  the  hill  were 
not  considered  for  travel  time  analysis. 

Figure  1.5  shows  the  parameters  and  ray  diagram  for  the  sea  floor  hill  model  which 
is  simply  a  layer  over  a  half-space.  For  each  super-critical  range  there  are  three  simple  paths 
for  compressional  energy  to  reach  the  OBH;  directly  through  the  water,  refracted  at  the 
critical  angle,  and  diffracted  from  the  hill.  These  account  for  the  first  three  arrivals  in  the 
data.  With  the  hill  acting  as  a  secondary  source,  simple  ray  tracing  and  Snell's  law  were 
used  to  create  the  travel  time  plots  of  figure  1.6.  Diffraction  curves  for  the  hill  at  various 
ranges  in  the  sagittal  plane  and  at  various  offsets  out  of  the  sagittal  plane  are  shown.  The 
DT s  are  very  similar  to  those  from  the  data  (figure  1.1).  The  diffracted  and  refracted 
phases  arrive  closest  in  time  at  a  range  of  approximately  7.2  km  for  a  hill  at  6  km.  in  figure 
1.6a.  At  points  two  kilometers  away  on  either  side  of  7.2  km.  (5.2  and  9.2  km.),  the 
offset  between  the  two  arrivals  has  increased  by  approximately  0.35  seconds  (DT).  The 
corresponding  DT  in  the  data  is  about  0.40  seconds. 

Since  the  position  of  the  hill  does  not  affect  the  direct  or  refracted  arrivals  (for  this 
simple  analysis)  the  traces  for  a  number  of  different  hill  locations  within  the  sagittal  plane 
are  shown  superimposed  in  one  plot  in  figure  1.6a.  As  expected,  the  shapes  of  the 
refraction  branch  diffractions  are  the  same  for  the  different  hill  locations.  Moving  the  hill 
away  from  the  sagittal  plane  (into  or  out  of  the  page  in  figure  1.5)  effects  mainly  the 
vertical  placement  of  the  diffraction  (in  time)  with  only  slight  changes  in  its  shape  (figure 
1 ,6b).  It  should  be  noted  that  the  model  geometry  used  in  these  simple  models,  as  well  as 
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in  the  buried  body  model  below,  is  the  same  as  the  data  geometry.  For  the  Finite  difference 
modeling,  the  two  geometries  are  not  the  same  and  care  must  be  taken  not  to  confuse 
them^. 


Ray  tracing  for  the  buried  body  model  was  done  using  velocity-depth  profile  A  of 
figure  1.4.  Again,  as  with  the  sea  floor  hill  model,  horizontal  movement  of  the  diffractor 
at  a  given  depth  within  the  sagittal  plane  does  not  alter  the  shape  of  the  diffraction,  only  the 
horizontal  placement  of  it  with  respect  to  the  first  refracted  arrival  (figure  1.7a).  However, 
the  shape  of  the  diffraction  can  be  changed  by  moving  the  diffractor  vertically  in  the  upper 
crust.  Moving  the  diffractor  deeper  tends  to  flatten  out  the  diffracted  arrival  shape,  as 
seen  in  figure  1.7b.  Vertical  movement  of  the  buried  body  within  the  upper  crust  also 
moves  the  diffracted  arrival  vertically  in  time.  If  the  body  is  much  below  the  sea  floor 
(more  than  about  300  meters),  the  diffraction  arrives  before  the  sea  floor  refraction  at  large 
ranges.  Of  course,  by  moving  the  body  away  from  the  sagittal  plane,  the  diffracted  arrival 
may  also  be  moved  vertically.  Because  of  this,  the  possible  depth  range  for  the  diffractor 
is  not  well  constrained.  However,  by  choosing  the  right  combination  of  diffractor  depth 
and  lateral  offset,  the  DTs  for  this  model  can  easily  be  made  to  match  those  seen  in  the 
data. 


Three  important  facts  have  been  brought  out  by  ray  tracing.  First,  and  most 
important,  is  that  assuming  diffractions  will  occur  from  a  hill  or  buried  body,  the  DTs  seen 
in  the  data  can  be  reproduced  by  the  simple  models  shown  above.  Minor  changes  in 
diffracted  arrival  shape,  from  depth  difference  or  different  velocity-depth  profiles,  do  not 
greatly  affect  the  DTs.  Second,  the  shape  of  the  diffraction  is  not  significantly  affected  by 
changes  in  horizontal  position(up  to  2  km)  of  the  diffracting  body  either  parallel  or 
perpendicular  to  the  sagittal  plane.  Because  of  this,  the  horizontal  placement  of  the 
diffractor  in  the  finite  difference  models  does  not  have  to  be  exactly  the  same  as  that  which 
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Figure  1.5.  Geometry  used  for  ray  tracing  models  (also  the  same  as  data  geometry). 
Three  ray  paths  from  each  shot  to  the  OBH  are  responsible  for  the  refracted,  diffracted 
and  direct  arrivals  seen  in  Figure  1.1. 
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Figure  1.6.  Travel  time  plots  created  from  ray  tracing  analysis  for  a  sea  floor  hill.  Figure 
a)  shows  h ills  at  different  ranges  in  the  sagittal  plane  of  the  shots  and  receiver.  Figure  b) 
shows  hills  at  a  fixed  range  (6  km)  and  offset  from  the  sagittal  plane  by  different  distances. 
The  shape  of  the  diffraction  is  not  significantly  affected  by  changes  in  the  position  of  the 
hill  (up  to  2  km)  in  either  case. 
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Figure  1.7.  Travel  time  plots  for  a  buried  body.  Figure  a)  shows  the  effect  of  moving  the 
body  horizontally  at  a  constant  depth  (0.4  km).  Figure  b)  shows  the  effect  of  moving  the 
body  vertically  at  a  given  range  (4  km). 
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produced  the  data.  Finally,  the  shape  of  the  diffraction  can  be  changed  by  altering  the 
velocity-depth  function  (the  shape  changes  between  the  hill  model  and  the  buried  body 
model  because  of  a  change  in  V-Z  function)  or  by  changing  the  depth  of  the  diffractor 
(which  effectively  changes  the  V-Z  function  for  the  scattered  energy). 
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FINITE  DIFFERENCE  MODELING 


Amplitude  and  phase  conversion,  as  well  as  travel  time,  of  the  seismic  energy  need  to 
be  reconstructed  in  order  to  get  a  complete  solution  of  the  problem.  This  section  discusses 
synthetic  seismograms  and  wavefront  'snapshots’  generated  by  a  finite  difference  solution 
of  the  elastic  wave  equation  for  heterogenous  media; 


pii  =  (X+p)V(V*u)  +  +  VX(V*u)  +  Vp.  x  (V  x  p)  +  2(Vp*V)u  (1) 


where  u  is  the  displacement  vector,  p  is  density,  X,p  are  Lame's  elastic  parameters,  V  is 
the  del  vector  operator  and  double  dots  signify  the  second  time  derivative.  Some  finite 
difference  schemes^  are  unstable  at  liquid-solid  boundaries  and  require  specifically  coded 
boundary  conditions.  A  differencing  scheme  developed  by  Nicole tis^  and  Bhasavanija^ 
is  used  for  the  solutions  in  this  study.  In  this  scheme,  the  elastic  (Lame's)  parameters  and 
density  at  a  point  are  calculated  from  defined  velocities  at  eight  surrounding  points.  In  this 
way,  a  smoothing  of  parameters  takes  place  at  boundaries  and  the  scheme  is  stable. 
Specific  boundary  conditions  for  the  liquid-solid  interface  are  therefore  not  needed  with 
this  formulation. 

Accuracy  and  stability  of  the  solution  depend  on  the  actual  differencing  scheme  used 
and  on  the  time  and/or  space  increments.  Stephen^  compared  his  finite  difference  solution 
to  the  reflectivity  method  of  Fuchs  and  Miiller^  and  found  that  when  boundary  conditions 
of  the  second  order  were  used  there  is  excellent  agreement  between  the  two  methods. 
Bhasavanijall  compared  his  results  to  Fourier  synthesis  waveforms  and  also  found 
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excellent  agreement.  If  the  space  increments  are  too  large,  the  seismic  energy  will  appear 
to  disperse,  with  low  frequencies  travelling  faster  than  high  frequencies^.  Kelly  et.  al 
found  that  dispersion  is  minimal  if  there  are  at  least  ten  grid  points  per  minimum  seismic 
wavelength.  Instabilities  caused  by  insufficient  time  sampling  can  be  avoided  if  the  time 
increment  satisfies; 

At  =  h  /  (0C“+p2)i/2  ,  oc  =  (tt.+24)/p)l/2  ,  (3=(p/p)i/2  (2) 


where  h  equals  the  space  increment  (horizontal  and  vertical  increments  are  equal  for  these 
models),  a  is  the  P-wave  velocity,  and  (3  is  the  S-wave  velocity. 

Figure  1.8  shows  the  set  up  of  the  general  finite  difference  grid.  Only  the  points  in  the 
transition  zone  are  calculated  using  the  heterogenous  finite  difference  scheme  of 
Bhasavanija.  The  layers  above  and  below  the  transition  zone  are  homogeneous  (both 
vertically  and  horizontally)  and  are  solved  using  a  homogeneous  solution.  The  source  is 
introduced  along  the  top  of  the  grid  and,  for  potential,  is  in  me  form  of  the  first  derivative 
of  a  Gaussian  curved.  Absorbing  boundaries'"^  are  used  at  maximum  range  and  depth 
and  there  is  an  axis  of  symmetry  at  zero  range. 

As  stated  earlier,  the  geometry  of  the  finite  difference  models  is  different  than  the 
geometry  of  the  data  and  ray  tracing  models.  For  the  data,  there  is  one  OBH  and  a  line  of 
shots  at  the  sea  surface  (figure  1.5).  The  situation  is  opposite  in  the  finite  difference 
models  with  one  surface  shot  and  a  line  of  receivers  at  the  sea  floor  (figure  1 ,9a)  (although 
the  vertical  position  of  the  receivers  can  be  varied).  Travel  time  plots  for  the  two  setups 
are  quite  different  for  a  hill  at  the  same  range  (compare  figures  1.6  and  1.9b).  Only  one 
trace  is  common  to  both  the  data  and  any  given  model  geometry,  depending  on  the  shot- 
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diffractor  separation.  To  actually  synthetically  reproduce  the  data,  one  model  geometry 
must  be  run  for  every  shot-receiver  arrangement  in  the  data.  By  changing  the  location  of 
the  diffractor  in  the  model,  a  different  shot-receiver  range  of  the  data  can  be  reproduced. 

A  Flat  Model 

A  control  model  with  a  flat  sea  floor  and  velocity-depth  profile  A  of  figure  1 .4  was 
first  run  to  produce  the  expected  arrivals  for  an  area  without  diffractors.  The  different 
phases  and  partitioning  of  energy  can  be  seen  in  wavefront  'snapshots’  (figure  1.10)  of 
divergence  and  curl  of  displacement.  Compressional  and  shear  energy  are  proportional  to 
the  squares  of  divergence  and  curl  of  displacement  respectively,  and  the  snapshots  can  be 
considered  as  representations  of  compressional  and  shear  wave  energy  travelling  across  the 
grid.  The  curl  amplitudes  have  been  scaled  by  (3)^/2  times  the  divergence  scaling  in  order 
to  best  represent  the  relationship  to  energy.  Snapshots  were  output  every  0.4  seconds  (320 
time  steps)  with  the  first  frame  at  two  seconds. 

All  of  the  phases  expected  for  the  FLAT  model  can  be  seen  in  figure  1.10.  The 
nomenclature  given  in  brackets  is  after  Brekhovskikh^  and  is  appropriate  for  the  phases 
present  in  the  case  of  a  sharp  interface  between  two  media  of  constant  velocities.  These 
symbols  are  not  exactly  correct  for  the  phases  seen  in  the  models  for  this  study  because  of 
the  presence  of  velocity  gradients  below  the  seafloor.  If  the  gradients  were  allowed  to 
decrease  to  zero,  the  phases  seen  in  the  snapshots  would  degenerate  to  those  given  in 
brackets.  Brekhovskikh’s  nomenclature  is  used  here  because  it  closely  approximates  the 
system  of  waves  which  exist  in  media  with  velocity  gradients  below  a  sharp  interface. 

The  source  is  located  2  km  above  the  top  of  the  grid  and  the  direct  water  wave  [Pi]  is 
just  above  the  seafloor  in  the  2.0  second  snapshot  (see  figure  1.10).  At  2.4  seconds  the 


45 


Figure  1.8.  Setup  of  the  grid  for  finite  difference  models.  Range  and  depth  increments 
are  equal  to  10  meters  (figure  not  to  scale).  Each  model  was  run  for  4000  time  steps  of 
0.00125  seconds  each.  The  source  is  introduced  along  the  top  boundary  of  the  grid.  The 
bottom  and  right  hand  side  boundaries  are  absorbing.  The  left  hand  boundary  is  an  axis  of 
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Figure  1.9.  Model  geometry  (figure  a)  and  associated  travel  time  plots  from  ray  tracing 
(figure  b) .  The  arrival  marked  virtual  hill  is  the  signal  from  the  opposite  side  of  the  axis 
of  symmetry  (axis  is  at  zero  range). 
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Figure  1.10a.  Wavefront  'snapshots'  of  divergence  and  curl  of  displacement  for  the 
FLAT  model.  Displacement  and  curl  represent  compressional  and  shear  energy 
respectively.  Amplitude  scaling  of  divergence  is  (3)^  times  that  of  the  curl  snapshots.  In 
addition,  the  range  of  amplitude  contoured  is  often  different  above  and  below  the  sea  floor 
and  is  given  on  the  right  side  of  each  frame  in  terms  of  the  amplitude  (in  percent  of 
maximum  amplitude)  which  is  contoured.  Only  the  direct,  reflected,  and  transmitted  waves 
are  seen  up  to  2.8  seconds.  (Nomenclature  in  brackets  after  Brekhovskikh^,  pp.  318- 
324). 
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direct  wave  has  reached  the  seafloor  and  partitioned  into  direct  [Pi  J,  transmitted  l P 1  P2l. 
and  reflected  [Pi Pi  ]  compressional  waves  and  a  transmitted  convened  shear  wave  f P 1 S 2) • 
If  the  lower  medium  were  a  uniform  half-space  the  transmitted  waves  (tPiP2l.  [ P 1 S 2 1) 
would  propagate  out  the  bottom  of  the  grid.  However,  in  the  presence  of  the  velocity 
gradients,  portions  of  these  waves  with  ray  parameter  greater  than  the  slownesses  at  the 
bottom  of  the  grid  (0.17  for  P  and  0.29  for  S)  are  refracted  in  the  upper  crustal  velocity 
gradients  and  impinge  on  the  seafloor  from  below.  These  are  referred  to  as  P  and  S  diving 
waves.  In  the  upper  layer  they  are  kinematically  similar  to  the  pure  head  waves  described 
in  Brekhovskikhi?  and  we  have  given  them  the  same  designation  ([P1P2P1].  [P1S2P1]). 
The  pure  head  waves  do  not  exist  in  our  model  because  they  are  only  defined  for  a 
homogeneous  lower  layer The  P-diving  wave  can  be  seen  in  the  water  column  as  it 
moves  out  in  front  of  the  direct  wave  after  3.2  seconds. 


When  the  P-diving  wave  is  incident  on  the  seafloor  from  below,  it  also  sets  up  a 
[P1P2S2]  converted  phase  (compressional  wave  through  the  water,  Pi,  P-diving  through 
the  upper  crust,  P2,  converted  shear  reflection  at  the  seafloor,  S2).  This  converted  S-wave 
is  coupled  to  the  P-diving  wave  at  the  seafloor,  and  is  also  tangential  to  the  transmitted  S- 
wave  at  depth.  At  3.6  seconds  the  S-diving  wave  [P1S2]  begins  to  move  away  from  the 
direct  wave  root  [(Pi ) l )  and  [(Pi )3].  Parentheses  indicate  a  wave  whose  amplitude 
decreases  exponentially  away  from  the  interface.  The  'direct  wave  root'  is  an  exponentially 
decaying  wave  in  the  lower  medium  corresponding  to  super-critical  reflections  in  the  upper 
medium^.  It  has  compressional  and  shear  parts,  [(Pi ) l ]  and  [(Pi )3]  respectively. 
Coupled  to  the  S-diving  wave  [P1S2]  is  a  P-wave  whose  amplitude  decreases 
exponentially  with  depth  [Pl(S2)]. 
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Interference  waves^,  which  are  essentially  the  superposition  of  diving  wave 
multiples  incident  on  the  seafloor  from  below,  occur  in  the  snapshots  at  3.2  and  4.0 
seconds  respectively.  After  they  have  separated  from  the  principle  diving  waves,  these 
waves  appear  very  much  like  pure  head  waves  in  the  water  column  and  are  referred  to  as 
interference  head  waves^.  The  interference  P-head  wave  can  be  seen  at  3.6  seconds  but 
the  interference  S-head  wave  cannot  be  seen  since  it  does  not  develop  as  a  distinct  phase 
until  after  the  last  snapshot  (later  than  4.4  seconds).  The  only  other  major  phase  to  occur 
is  the  pseudo-Rayleigh  wave  which  appears  in  the  4.4  second  snapshot.  This  phase  is  an 
evanescent  wave  in  the  crust,  with  both  divergence  and  curl  componenets,  and  it  appears 
as  a  plane  propagating  wave  in  the  water. 

Along  with  these  expected  arrivals,  there  are  a  few  phases  in  the  snapshots  which  are 
artifacts  of  the  grid  and  the  numerical  method.  The  P-diving  wave  reflects  from  the  top  of 
the  grid  (seen  in  the  4.0  second  snapshot  of  figure  1.10)  because  this  boundary  is  rigid 
after  the  source  is  introduced.  Also,  there  is  some  noise  in  the  later  snapshots  along  the 
bottom  of  the  grid  due  to  incomplete  absorption  of  energy  along  this  boundary.  The 
phases  in  the  water  column  of  the  curl  snapshots  are  artifacts  of  reflection  from  the  upper 
rigid  boundary.  Water  column  wavefronts  which  occur  in  the  snapshots  of  the  curl  of 
displacement  (figure  1.10 )  would  not  be  present  if  shear  energy  was  plotted.  Since  the 
shear  velocity  in  water  is  zero,  the  amplitude  of  the  shear  energy  in  the  water  would  also 
be  identically  zero.  P-S  conversion  below  the  water-solid  interface  is  accurate  and 
complete.  None  of  these  problems  interfere  with  any  phases  of  interest  for  this  study. 

Figure  1.11  shows  the  time  series  of  pressure  for  receivers  at  depths  of  2.48  and  2.98 
km.  and  of  vertical  displacement  for  receivers  at  3.48  km  depth.  The  time  series  in  figures 
1 .11a  and  1 ,11b  are  very  similar  except  for  slight  offsets  of  the  arrivals  in  time  and  space 
(between  the  two  seismograms)  and  the  number  of  water  wave  multiples  The  seafloor 
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Figure  1.10b.  Continuation  of  the  Wavefront  'snapshots'  of  divergence  and  curl  of 
displacement  for  the  FLAT  model  seen  in  figure  1.10a.  Note  that  the  range  window 
moves  with  the  wavefronts.  At  3.2  seconds  the  P-diving  wave  can  be  seen  as  a  plane 
wave  in  the  water  and  the  direct  wave  root  ^  travels  below  the  sea  floor  directly  beneath 
the  direct  wave.  (Nomenclature  in  brackets  after  Brekhovskikh^,  pp.  318-324). 
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Figure  1.10c.  Continuation  of  the  wavefront  'snapshots'  of  divergence  and  curl  of 
displacement  for  the  FLAT  model  seen  in  figures  10a  and  10b.  Note  that  the  range 
window  moves  with  the  wavefronts.  The  planar  waves  in  the  water  column  which 
correspond  to  the  pseudo-  Rayleigh,  S-diving,  and  interference  P-head  waves  can  be  seen 
as  separate  phases  in  the  4.4  second  snapshot.  (Nomenclature  in  brackets  after 
Brekhovskikh  17,  pp.  318-324). 
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Figure  1.11.  Synthetic  seismograms  for  the  FLAT  model.  Figure  a)  is  the  time  series  of 
pressure  for  receivers  in  the  water  column  (a  depth  of  2480  m  or  520  meters  above  the 
seafloor).  Figure  b)  is  the  time  series  of  pressure  for  receivers  near  the  ocean  bottom  (20 
meters  above  the  seafloor).  Figure  c)  is  the  time  series  of  vertical  displacement  for 
receivers  in  the  bottom  (a  depth  of  3480  m  or  480  m  below  the  sea  floor). 
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receivers  (figure  1.11b)  record  only  water  waves  (direct  and  multiples)  incident  from  above 
but  the  receivers  in  the  water  column  (figure  1.11a)  record  both  up  and  downgoing  waves 
(direct  and  multiples).  Therefore,  figure  11a  has  twice  the  number  of  multiples  of  1  .lib.  As 
in  the  snapshots,  these  multiples  are  reflections  from  the  top  of  the  grid  rather  than  from 
the  true  sea  surface. 

Ranges  and  times  of  the  arrivals  seen  in  the  seismograms  correspond  to  those  seen  in 
the  snapshots  of  figure  1.10.  Major  arrivals  seen  in  the  water  column  (figure  1.11a)  and 
seafloor  (figure  1.11b)  seismograms  are  the  P-diving  and  interference  head  waves,  the 
converted  S-diving  wave  and  the  pseudo-  Rayleigh  wave.  The  pseudo-Rayleigh  wave 
looks  like  another  plane  wave  in  the  water. 

A  few  other  arrivals  appear  in  figure  1.11c,  the  vertical  displacement  seismograms  for 
buried  receivers  (480  m  below  the  seafloor).  The  P  and  S-diving  waves  appear  as  they  do 
in  the  other  seismograms,  and  the  pseudo-  Rayleigh  wave  is  present  but  of  a  low 
amplitude  since  its  amplitude  decays  exponentially  with  depth.  Interference  P  and  S  waves 
immediately  follow  the  P  and  S-diving  waves  respectively.  The  P1P2S2  arrival,  which  has 
no  corresponding  water  column  phase,  can  be  seen  in  the  seismogram  arriving  after  the 
interference  P-wave  (at  ranges  less  than  about  8  km)  and  before  the  S-diving  wave.  It 
becomes  indistinguishable  from  the  interference  P  wave  at  large  ranges.  Figure  1.11c  also 
contains  arrivals  from  the  direct  wave  root,  which  is  the  crustal  expression  of  the  direct 
water  wave  seen  in  figure  1.11b. 

Diffraction  Models 

Five  models  with  diffractors  were  run  to  determine  if  diffractions  would  occur,  and  if 
so,  to  look  at  the  effect  that  size  and  shape  of  the  diffracting  body  has  on  the  diffractions. 
The  anomalous  bodies  (figure  1.12)  were  placed  at  a  range  of  five  kilometers.  This  gives 
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Figure  1.12.  Relative  sizes  and  shapes  of  different  diffractors  used  for  the  various  models. 
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:he  seismic  energy  sufficient  time  to  partition  in  order  to  examine  the  effects  of  different 
phases  hitting  the  body.  Four  models,  HILL,  VALLEY1,  VALLEY2,  and  BLOCK,  had 
relatively  small  diffractors  and  used  velocity-depth  profile  A  (figure  1.4).  VALLEY3,  the 
final  model  run,  had  a  much  larger  diffractor  and  used  velocity-depth  profile  B. 

Models  HILL,  VALLEY1,  VALLEY2,  and  BLOCK  ail  produced  definite  diffractions, 
an  example  of  which  can  be  seen  in  figures  1.13  and  1.14  for  model  VALLEY1.  The 
snapshots  for  VALLEY1  (figure  1.13)  start  at  3.6  seconds.  The  earlier  snapshots  are  not 
presented  because  of  their  similarity  to  the  early  snapshots  of  the  FLAT  model  (figure 
1.10,  which  can  be  used  for  reference).  Three  features  distinguish  VALLEY1  from  FLAT 
and  can  be  seen  in  the  corresponding  snapshots  (figure  1.10  for  FLAT  and  figure  1.13  for 
VALLE Yl).  First,  the  P-diving  wave  in  the  water  column  is  distorted  by  the  presence  of 
the  valley.  Second,  when  the  P-diving  wave  hits  the  valley,  it  produces  a  diffracted  S- 
diving  wave  in  the  crust  (seen  best  in  the  4.0  and  4.4  second  snapshots).  A  P-diffracted 
wave  probably  also  occurs  but  is  indistinguishable  from  the  interference  head  wave. 
Finally,  and  most  important,  when  the  direct  wave  and  its  root  hit  the  valley,  definite 
diffracted  compressional  and  shear  waves  occur  and  radiate  out  from  the  valley  (seen  in  the 
4.0  and  4.4  second  snapshots).  The  seafloor  seismogram  for  VALLEY1  (figure  1.14) 
shows  all  of  the  arrivals  seen  in  the  FLAT  model  (figure  1.1 1)  as  well  as  the  S-difffacted 
arrival  (from  P-diving)  and  the  P-diffraction  (from  the  direct  wave).  The  S-diffraction  set 
up  by  the  direct  wave  is  not  distinguishable  from  the  pseudo-Rayleigh  wave  in  the 
seismogram.  Although  the  snapshots  and  seismograms  for  models  BLOCK,  HILL,  and 
VALLEY2  are  not  presented  here,  the  same  three  effects  are  present  in  all  of  these  models. 

VALLEY3  failed  to  produce  a  significant  diffraction  (see  figure  1.15).  The 
expected  P-diving  and  interference  head  waves  and  the  S-diving  wave  are  all  seen  in  this 
model  as  in  the  FLAT  model.  The  P-diving  wave  is  again  distorted  by  the  presence  of  the 
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valley  but  no  significant  diffractions  occur  either  when  the  P-diving  or  the  direct  wave  hit 
the  valley.  Because  the  upper  crustal  velocities  are  lower  in  this  example  (V-Z  function  B) 
than  in  the  previous  examples  (V-Z  function  A),  the  P-diving  and  interference  waves  are 
much  more  separated. 

Figure  1.16  shows  the  relative  amplitudes  of  the  diffracted  arrivals  and  other  phases 
for  the  seafloor  trace  at  7  kilometers  range  of  all  six  models.  Using  the  trace  from  the 
FLAT  model  for  reference,  it  is  obvious  that  the  shape  of  the  diffractor  is  very  important  in 
determining  the  character  of  the  signals.  VALLE Y1  and  VALLE Y2  have  almost  identical 
signals  except  for  slight  differences  in  amplitudes.  The  BLOCK  and  HILL  traces  are 
different  from  each  other  and  from  the  V  ALLEY  1  and  VALLEY2  models.  The  trace  from 
VALLEY3  is  different  from  all  of  the  others  not  only  because  of  the  size  of  the  valley  but 
also  because  >t  was  produced  with  a  different  velocity -depth  function.  All  of  the  models 
which  contain  significant  diffractions  have  the  same  arrivals,  but  often  with  different 
arrival  times.  The  amplitudes  of  the  diffractions,  as  well  as  most  of  the  other  arrivals,  are 
of  the  same  magnitude  from  one  model  to  the  next.  The  only  phase  which  varies 
appreciably  between  the  models  is  the  interference  head  wave  which  appears  to  be 
significantly  affected  by  topography. 
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Figure  1.13.  Divergence  and  curl  of  displacement  snapshots  for  the  VALLEY1  model 
(corresponding  to  compressional  and  shear  energy).  Earlier  timesteps  are  very  similar  to 
those  of  the  FLAT  model.  Distortion  of  the  P-diving  wave  can  be  seen  at  3.6  seconds. 
The  diffracted  phases  from  the  P-diving  wave  and  the  direct  wave  begin  to  appear  in  the 
4.0  second  snapshot  and  are  quite  clear  as  separate  phases  at  4.4  seconds. 
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Figure  1.14.  Synthetic  seismogram  for  seafloor  receivers  in  the  VALLEY1  model.  The 
seismogram  is  very  similar  to  that  of  the  FLAT  model  (see  figure  1.11)  except  for  the  two 
diffracted  arrivals  and  the  amplitude  scaling.  A  greater  amplitude  scaling  was  used  here  to 
accentuate  the  diffracted  arrival.  Other  arrivals  in  the  seismogram  have  approximately  the 
same  amplitude  as  in  figure  1.11.  The  P-diffraction  from  the  direct  wave  intercepts  the 
direct  wave  branch  and  has  a  velocity  (slope)  similar  to  the  P-diving  wave.  The  P-S 
diffraction  is  tangential  to  the  P-diving  wave  and  has  a  velocity  similar  to  the  S-diving 
wave. 
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Figure  1.15a.  Divergence  and  curl  of  displacement  snapshots  (corresponding  to 
compressional  and  shear  energy)  for  ihe  model  VALLEY3.  Amplitude  scaling  of  these 
snapshots  is  0.5  that  of  figures  1.10  and  1.13.  Note  also  that  the  range  window  moves 
with  the  wavefronts.  The  water  column  plane  wave  corresponding  to  the  P-  diving  wave 
is  distorted  by  the  presence  of  the  valley  in  the  3.6  and  4.0  second  snapshots.  However, 
diffracted  phases  are  indiscernible  in  any  of  the  snapshots.  The  presence  of  prominent  S- 
wave  phases  in  the  curl  snapshots  indicates  that  strong  P-S  conversion  takes  place  with  the 
velocity-depth  function  of  this  model  (profile  B  of  figure  1.4). 
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Figure  1.15b.  Continuation  of  the  divergence  and  curl  of  displacement  snapshots 
(corresponding  to  compressional  and  shear  energy)  for  the  model  VALLEY3  as  seen  in 
figure  1.15a.  Amplitude  scaling  of  these  snapshots  is  0.5  that  of  figures  1.10  and  1.13. 
Note  also  that  the  range  window  moves  with  the  wavefronts.  The  water  column  plane 
wave  corresponding  to  the  P-  diving  wave  is  distorted  by  the  presence  of  the  valley  in  the 
4.4  second  snapshot.  However,  diffracted  phases  are  indiscernible  in  any  of  the 
snapshots.  The  presence  of  prominent  S-wave  phases  in  the  curl  snapshots  indicates  that 
strong  P-S  conversion  takes  place  with  the  velocity-depth  function  of  this  model  (profile  B 
of  figure  1.4). 
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Figure  1. 16.  Comparison  of  the  7  km  trace  for  seafloor  receivers  in  each  of  the  six  models. 
Amplitude  gain  for  the  top  traces  is  ten  times  that  of  the  bottom  traces.  The  P-diving  (P), 
S-diving  (S),  and  pseudo-Rayleigh  (pR)  wave  arrivals  are  similar  for  all  of  the  models  (the 
delay  in  arrival  times  for  model  VALLEY3  is  caused  by  a  different  velocity-depth  profile). 
The  relative  arrival  times  of  the  diffracted  P  from  the  direct  wave  root  (dP)  and  the 
diffracted  S  from  the  P-diving  wave  (dS) ,  and  the  amplitude  of  the  interference  P-head 
wave  (IHW)  change  between  the  models.  Note  the  lack  of  diffracted  arrivals  for  model 
VALLEY3. 
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Phase  identification  was  accomplished  by  a  number  of  relatively  simple  techniques. 
Snapshots  of  divergence  and  curl  of  displacement  (proportional  to  compressional  and 
shear  energy  respectively)  are  quite  effective  at  showing  the  partitioning  of  energy  between 
P  and  S  waves.  Another  way  to  identify  phases  in  the  snapshots  or  time  series  is  to  plot 
their  particle  motion.  Figure  1.17  shows  particle  motion  plots  for  the  major  phases 
present  in  the  time  series  plots  for  the  FLAT  model  (figures  1.10  and  1.11).  These  plots 
confirm  that  panicle  motion  for  P-waves  is  along  the  direction  of  wave  propagation.  S- 
wave  motion  is  transverse  to  propagation  direction  and  pseudo-Rayleigh  wave  motion  is 
retrograde  elliptical  (at  the  seafloor). 

The  interference  wave  is  essentially  the  result  of  the  superposition  of  all  possible 
multiples1^.  When  this  wave  has  separated  from  the  principle  diving  wave  it  appears 
similar  to  a  pure  head  wave  in  the  water  column  and  is  referred  to  as  an  interference  head 
wave.  Both  P  and  S  interference  waves  appear  in  the  models  shown  here  (see  figure  1.10) 
but  only  the  interference  P-head  wave  totally  separates  from  the  principle  P-diving  wave  in 
the  time  range  available.  The  horizontal  phase  velocity  of  the  interference  P-head  wave  is 
between  the  velocities  of  the  P-diving  wave  and  the  pure  P-head  wave  (which  exists  for  a 
half-space  with  uppermost  P-velocity).  In  this  model,  for  the  ranges  present,  the  velocity 
of  the  interference  P-head  wave  is  dominated  by  the  first  and  second  multiple.  From  ray 
theory,  these  should  separate  from  the  interference  head  wave  at  ranges  of  approximately 
6.5  and  8.5  km  respectively  but  are  of  very  low  amplitude  and  are  indistinguishable  from 
the  interference  head  wave.  Cerveny  and  Ravindra^  predict  that  the  diving  wave  and  the 
interference  head  wave  will  dominate  the  seismogram  and  the  separated,  non-interference 
multiple  diving  waves  will  be  weaker.  Our  results  support  the  predictions  of  Cerveny  and 
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Ravindra.  However,  the  velocity  gradients  used  for  this  study  (4.0  sec.'l  at  the  seafloor) 
are  much  greater  than  the  those  used  by  Cerveny  and  Ravindra  (<  0.04  sec.’l).  Data  from 
OBH  2W  lack  any  interference  waves,  because  topography  along  the  line  causes 
attenuation  due  to  scattering.  Indeed,  the  amplitude  and  coherence  of  the  interference  P- 
head  wave  arrival  in  the  HILL  model  is  decreased  from  the  amplitude  in  the  FLAT  model 
simply  by  adding  a  single  hill  onto  the  sea  floor  (compare  figures  1.10  and  1.18). 

Stephen^  discussed  the  generation  of  'double  head  waves'  due  to  interaction  with  a 
seafloor  hill.  In  this  example  the  compressional  wave  critical  angle  is  reached  twice  by  a 
direct  wave  travelling  over  the  hill  (once  on  the  slope  of  the  hill  and  again  on  the  flat 
surface  after  it).  The  situation  is  somewhat  more  complicated  when  a  velocity  gradient  is 
present  since  diving  waves  rather  than  head  waves  occur  and  there  is  no  single  critical 
angle  for  diving  waves.  Double  head  waves  do  not,  however,  occur  in  the  models  of  this 
study,  because  the  anomalous  structures  were  placed  far  beyond  the  critical  ranges  for  any 
diving  compressional  energy.  However,  when  the  P-diving  wave  hits  the  diffractor,  it 
causes  a  convened  P-S  diffracted  phase  to  occur.  Since  this  phase  occurs  in  the  BLOCK, 
a.s  well  as  the  HILL  and  VALLEY  models,  it  appears  that  it  must  be  a  diffraction  rather 
than  some  son  of  ’double  head  wave'  refraction. 

The  occurrence  of  backscattered  energy  also  verifies  the  hypothesis  that  diffraction 
r.  her  than  refraction  or  reflection  of  seismic  energy  is  occuring  around  the  seafloor 
fe  ttures.  Backscattering  of  energy  is  most  convincingly  seen  in  figure  1.19,  the  time 
se  les  of  pressure  for  a  buried  receiver  in  the  BLOCK  model.  Arrivals  predicted  by  the  flat 
model,  as  well  as  the  P-  diffraction  (from  the  direct  wave)  are  all  present  in  the 
seismogram.  Since  this  is  the  pressure  time  series,  S-waves  do  not  have  large  amplitudes 
(compare  this  with  the  vertical  displacement  seismogram  for  buried  receivers  of  the  FLAT 
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Figure  1.17.  Panicle  motion  plots  for  major  arrivals  of  the  FLAT  model.  The  large  arrow 
indicates  the  direction  of  wave  propagation  and  the  cross  indicates  the  onset  of  motion  for 
each  plot.  Receivers  for  1.17a-e  are  buried  at  480  m  below  the  seafloor  and  for  I.17f  the 
receivers  are  at  the  seafloor.  The  P-diving  wave  (1.17a)  and  the  interference  P-wave 
( 1 . 17b)  both  have  particle  motion  parallel  to  propagation  direction.  The  S-diving  wave 
(1.17d),  P1P2S2  (1.17c),  and  the  interference  S-wave  (1.17e)  all  have  panicle  motion 

perpendicular  to  propagation  direction  and  the  pseudo- Rayleigh  wave  (1.170  has 
retrograde  elliptical  panicle  motion  at  the  seafloor.  All  of  the  plots  are  for  receivers  at  a 
range  of  7  km. 
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Figure  1.18.  Divergence  of  displacement  snapshot  from  model  HILL  showing  almost 
complete  lack  of  interference  P-head  wave  (compare  this  with  4.4  sec  snapshots  of  Figures 
1.10  and  1.13). 
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Figure  1.19.  Synthetic  seismograms  for  buried  receivers  in  the  BLOCK  model  showing 
backscattering  of  energy  from  the  diffractor. 
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model,  figure  1.11c).  Only  the  pressure  from  P-waves  coupled  to  the  S- waves  will  show 
up  in  the  pressure  plot.  The  S-diffraction  from  the  P-  diving  wave  is  below  the  noise 
level  of  the  seismogram  for  this  reason.  The  signal  of  the  diffraction  is  quite  similar  to  that 
predicted  by  ray  tracing  for  a  buried  body  using  the  modf 1  geometry  (figure  1 .9b).  The 
two  limbs  of  the  diffraction  are  symmetric  about  the  range  of  the  diffractor  (5  km  for  the 
BLOCK  model)  as  in  figure  1.9. 

The  relationship  between  seismic  wavelength  and  size  of  the  diffractor  is  complicated 
by  a  number  of  factors.  Since  the  diffrators  modeled  in  this  study  are  at  the  water-crust 
interface,  it  is  not  immediately  obvious  if  the  direct  or  transmitted  energy  is  causing  the 
diffraction.  Velocity  gradients  cause  the  wavelength  of  the  seismic  energy  to  change 
rapidly  with  depth  and  the  wavefront  will  interact  differently  with  different  sizes  of 
diffractors  depending  on  depth.  It  appears  that  the  diffractions  are  the  result  of  subsurface 
energy.  The  P-diving  wave  is  incident  from  below  when  it  hits  the  diffractors  and  the 
subsurface  energy  for  the  direct  wave  is  in  the  form  of  its  prominent  direct  wave  root. 

Also,  a  diffraction  occurs  in  the  BLOCK  model  which  has  no  surface  topography  (no  hill 
or  valley)  to  influence  the  energy  actually  travelling  in  the  water.  If  it  is  subsurface  energy 
causing  the  diffractions,  then  the  sizes  of  the  diffractors  must  be  similar  to  the  P- 
wavelength  in  the  upper  crust.  Indeed,  if  the  wavefront  sees  only  one  half  of  the  diffractor 
at  a  time  (i.e.  it  first  encounters  the  upslope  of  a  hill  anad  later  encounters  the  downslope), 
then  the  half  width  of  the  diffractors  is  almost  identical  to  the  average  P-wavelength  of  the 
upper  crust  (half  width  of  450  meters  for  HILL,  VALLEY1  and  BLOCK  models  and  an 
average  P-wavelength  of  400  meters  for  the  upper  400  meters  of  crust).  For  model 
VALLEY3,  which  has  a  diffractor  of  half  width  equal  to  1320  meters  and  an  average  P- 
wavelength  of  320  meters,  no  diffraction  occurs  apparently  because  the  half  width  is  much 
larger  than  the  seismic  wavelength. 
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The  possibility  of  diffractions  occuring  from  a  number  of  different  anomalous  bodies 
has  been  clearly  established.  On  one  hand,  the  smaller  diffractors  of  HILL,  BLOCK, 
VALLEY1,  and  VALLEY2  caused  diffractions  but  are  generally  not  seen  in  the 
corresponding  bathymetry  (figure  1.4).  On  the  other  hand,  VALLEY3,  a  closer 
representation  of  the  bathymetry  of  the  data  area,  does  not  cause  a  diffraction.  One 
possible  compromise  to  the  problem  is  that  the  hill  offset  behind  the  valley  in  figure  1.4 
could  be  the  source  of  the  diffracted  energy.  Because  it  is  offset  behind  a  large  valley  its 
relative  size  is  closer  to  the  size  of  the  diffractors  in  the  earlier  models.  The  other  hills  in 
the  section  are  too  low  and  broad  to  act  as  diffractors. 

Assuming  that  the  offset  hill  can  be  the  source  of  the  diffracted  energy,  a  number  of 
comparisons  between  the  data  and  the  models  can  be  made.  If  the  difference  in  geometry' 
between  model  and  data  is  considered,  the  travel  times  for  the  diffraction  in  the  models 
(HILL,  VALLEY1,  VALLEY2,  and  BLOCK)  are  approximately  the  same  as  those  in  the 
data.  The  amplitudes  of  the  diffracted  arrivals,  while  not  as  large  as  the  first  P-refraction 
arrival,  are  significant.  It  is  interesting  that  in  the  data,  the  diffraction  is  even  larger  in 
amplitude  fhan  the  P-diving  wave  for  some  ranges.  Presumably,  both  the  travel  times  and 
amplitudes  of  the  model  diffractions  could  be  altered  by  perturbation  of  the  model 
velocity-depth  functions. 

Dirtct  comparison  of  the  data  and  models  is  not  possible  because  of  differences  in 
geometry  and  source  signatures  for  the  two  setups.  We  stress  that  the  objective  of  this 
study  was  to  demonstrate  the  phenomenon  of  re  fraction  branch  diffraction,  and  not  to 
match  exact  waveforms  in  the  data.  Simple  ray  tracing  showed  that  movement  of  the 
difffactor  did  not  significantly  change  the  shape  of  the  diffracted  arrival.  We  chose  to  place 
the  diffractors  at  5  km  range  in  the  models  to  best  take  advantage  of  the  finite  difference 
grid  without  regard  to  the  exact  matching  of  any  one  data  trace  geometry’. 
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The  signatures  of  the  synthetic  and  airgun  sources  are  different  enough  to  make 
waveform  matching  between  the  models  and  data  unreasonable.  The  synthetic  source  is  a 
single  pulse  with  peak  frequency  at  10  hz  (see  figure  1.20).  Airgun  sources  typically  have 
an  initial  pressure  pulse  plus  a  number  of  secondary  bubble  pulse  reverberations.  A  single 
pulse  source  was  used  for  the  modeling  in  order  to  better  define  and  identify  the  diffracted 
arrival.  The  exact  arrival  of  the  refraction  branch  diffraction  in  any  one  data  trace  is  often 
masked  by  the  airgun  source  signature  of  the  refraction  branch  of  the  seismogram  (figure 
1.1).  Only  when  the  entire  survey  of  OBH  2W  is  taken  as  a  whole  does  the  refraction 
branch  diffraction  stand  out  as  a  coherent  arrival.  The  synthetic  source  used  for  the  finite 
difference  modeling  allows  the  diffracted  arrival  to  be  picked  much  more  easily  on  any 
single  trace  (figure  1.16). 

A  complicating  factor  which  occurs  in  the  models  but  not  in  the  data  is  the  appearance 
of  shear  and  pseudo-Rayleigh  waves.  Direct  wave  energy  lost  to  shear  energy  by  P-S 
conversion  couiu  potentially  be  a  source  of  the  additional  energy  needed  to  increase  the 
diffraction  amplitudes.  White  and  Stephen^,  showed  that  P-3  conversion  is  decreased  by 
increasing  the  Poisson’s  ratio  of  the  basement  or  by  increasing  the  impedance  contrast 
between  the  S-wave  in  basement  and  the  P-wave  velocity  in  the  overlying  layer.  Spudich 
and  Helmberger^l  produced  two  ocean  bottom  models  which  do  not  contain  shear  waves, 
one  which  has  a  crust  of  high  Poisson's  ratio  (ratio  of  0.38  rather  than  0.25  for  a  Poisson 
solid)  and  the  other  which  had  a  high  impedance  contrast  at  the  top  of  the  basement  (crustal 
S-wave  velocity  of  2.023  km/sec  and  sediment  P-wave  velocity  of  3.0  km/sec).  The 
young  crust  in  the  ROSE  area,  however,  is  unsedimented  and  has  water  directly  on  top  of 
basement.  It  is  not  clear  how  the  substitution  of  a  water  layer  above  the  basement  for  the 
Spudich  and  Helmberger  transition  zone  or  White  and  Stephen  sediment  layer  would  affect 
P-S  conversion  at  the  interface. 
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From  the  me  'elir.g  done  in  this  study  it  is  dear,  however,  that  simply  lowering  the 
upper  crustal  P-wave  velocities,  while  assuming  a  Poisson  solid,  does  not  eliminate  P-S 
conversion  (conversion  occurred  for  both  velocity-depth  functions  of  figure  1.3).  A  more 
realistic  model  of  crustal  Poisson's  ratios  (rather  than  assuming  a  uniform  Poisson's  ratio 
of  0.25)  would  be  one  with  some  depth  dependence  caused  by  greater  fracturing  and 
porosity  of  the  crust  at  shallow  depths.  A  Poisson's  ratio  closer  to  0.35  would  be  more 
realistic  for  young,  unsedimented  crust  found  in  this  part  of  the  ROSE  area.  Other  factors 
which  affect  the  appearance  of  S -waves  in  marine  seismic  data  are;  angle  of  incidence  of 
the  P-wave,  incoherency  of  S- arrivals  due  to  rough  topography  and  lateral  crustal 
heterogeneities,  and  absorptioi  of  shear  energy  (higher  absorption  in  fractured,  porous 
crust). 


The  strength  of  P-S  conversion,  while  it  does  net  fit  with  the  data,  is  fortuitous  in 
that  it  demonstrates  the  pseudo-Rayleigh  wave  at  large  times  and  ranges  (see  figures  1.10 
and  1.17).  There  is  a  definite  break  between  the  shear  and  pseudo-Rayleigh  waves  and  the 
amplitudes  of  the  two  are  comparable  (figure  1.16).  Ray  parameter  analysis  on  the 
wavefront  shows  that  it  is  travelling  with  a  velocity  of  2.0  km/sec.  If  the  pseudo-Rayleigh 
wave  travels  at  92%  of  the  shear  wave  velocity  (the  free  surface  Rayleigh  wave  velocity 
relationslrp1)  then  this  wave  corresponds  to  an  S-wave  of  velocity  2.17  km/sec  at  an 
effective  depth  of  145  meters  instead  of  being  affected  by  the  interface  S-wave  velocity  of 
1.85  km/sec.  Rather  than  travelling  at  the  velocity  right  at  the  interface,  the  velocity 
grad;,\.  cauoes  the  pseudo-Rayleigh  phase  to  travel  at  an  'effective'  velocity  which  is 
related  to  the  seismic  wavelength  and  the  shear  velocity  gradient. 
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Figure  1.20.  Normalized  pressure  signatures  of  the  synthetic  and  airgun  array  sources. 

The  synthetic  source  is  the  third  derivative  of  a  Gaussian  curve  and  has  a  peak  frequency 
of  10  hz  with  an  upper  half  power  frequency  of  13.5  hz22.  The  seafi0or  reflection  pulse  is 
from  the  Lamont  Doherty  Geophysical  Observatory  airgun  array  used  for  ROSE  airgun  line 
1  (J.  Diebold,  pers.  comm.). 
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CONCLUSIONS 


The  following  conclusions  can  be  made  from  this  study: 

1 .  A  diffraction  has  been  observed  in  data  from  one  OBH  of  the  Rivera  Ocean  Seismic 
Experiment.  We  call  this  arrival  the  "refraction  branch  diffraction"  because  of  its  location 
on  the  refraction  branch  at  large  ranges  on  the  seismogram,  and  to  distinguish  it  from 
diffractions  which  occur  at  near  normal  incidence. 

2.  Finite  difference  modeling  of  seafloor  heterogeneities  has  shown  that  diffractions  w  ill 
occur  around  hills  and  valleys  of  the  size  commonly  found  on  the  ocean  bottom.  It  is  also 
possible  for  subsurface  bodies  of  large  velocity  contrast  to  cause  diffractions.  The  source 
of  the  diffractions  is  subsurface  energy  in  the  form  of  the  direct  wave  root  or  P-diving 
waves.  The  half-width  of  the  diffracting  bodies  are  all  similar  to  the  P-wavelength  in  the 
upper  crust. 

3.  The  "refraction  branch  diffraction"  in  the  data  from  OBH  2W  probably  originates  from 
the  hill  slightly  offset  from  the  track  line.  The  valley  seen  at  large  range  is  too  broad  to 
cause  a  diffraction  but  its  low  relief  would  accentuate  the  size  of  the  hill  behind  it,  allowing 
the  hill  to  act  as  a  distinct  diffractor.  Because  of  the  occurence  of  this  hill,  a  buried  block  is 
not  required  to  explain  the  diffraction.  Other  diffractions  may  occur  along  airgin  line  1  but 
destructive  interference  caused  by  rough  topography  may  ruin  their  coherence. 

4.  P-S  conversion  at  the  seafloor  is  prevalent  in  the  finite  difference  models  run  for  this 
study  but  not  in  the  data  from  OBH  2W.  Using  a  low  upper  crustal  P-wave  velocity  and 
assuming  a  Poisson's  solid  does  not  eliminate  the  shear  conversion.  If  the  velocity-depth 
function  of  Ewing  and  Purdy^  (V-Z  function  B  in  figure  1.3)  is  reasonable  for  the  upper 
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crust  then  it  is  apparent  that  the  assumption  of  a  Poisson's  ratio  of  0.25  is  not  valid.  A 
Poisson's  ratio  closer  to  0.35  is  proposed  as  more  realistic  for  this  region  but  more 
analysis  needs  to  be  done  to  correctly  determine  the  depth  dependence  of  Poisson's  ratio. 

5.  A  pseudo-Rayleigh  wave  occurs  at  large  ranges  using  the  V-Z  functions  appropriate  for 
the  ROSE  area  and  assuming  a  Poisson  solid.  The  wave  travels  at  the  interface  at  92%  of 
the  S-wave  velocity  at  an  'effective'  depth  related  to  the  seimic  wavelength  and  the  shear 
velocity  gradient  used  (approximately  145  meters  for  the  models  with  V-Z  function  A). 

6.  Interference  waves,  seen  in  the  models  but  not  in  the  data  from  OBH  2W,  are 
attenuated  by  scattering  from  rough  topography  along  the  line  and  are  below  the  noise 
level  in  the  seismograms. 
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Seismic  energy  partitioning  and  scattering  in  laterally 
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Seismic  energy  partitioning  and  scattering  in  laterally 
heterogeneous  ocean  crust 


ABSTRACT 

We  present  finite  difference  forward  models  of  elastic  wave  propagation  through 
laterally  heterogeneous  upper  oceanic  crust.  The  finite  difference  formulation  is  a  2-D 
solution  to  the  elastic  wave  equation  for  heterogeneous  media  and  implicitly  calculates  P 
and  SV  propagation,  compressional  to  shear  conversion,  interference  effects  and  interface 
phenomena.  Random  velocity  perturbations  with  Gaussian  and  self-similar  autocorrelation 
functions  and  different  correlation  lengths  (a)  are  presented  which  show  different 
characteristics  of  secondary  scattering.  Heterogeneities  scatter  primary  energy  into 
secondary  body  waves  and  secondary  Stoneley  waves  along  the  water-solid  interface.  The 
presence  of  a  water-solid  interface  in  the  models  allows  for  the  existence  of  secondary 
Stoneley  waves  which  account  for  much  of  the  seafloor  'noise'  seen  in  the  synthetic 
seismograms  for  the  laterally  heterogeneous  models. 

'Random'  incoherent  secondary  scattering  generally  increases  as  ka  (wavenumber, 
k,  and  correlation  length,  a)  approaches  one.  Deterministic  secondary  scattering  from 
larger  heterogeneities  is  the  dominant  effect  in  the  models  as  ka  increases  above  one. 
Secondary  scattering  also  shows  up  as  incoherence  in  the  primary  traces  of  the 
seismograms  when  compared  to  the  laterally  homogeneous  case.  Cross-correlation 
analysis  of  the  initial  P-diving  wave  arrival  shows  that,  in  general,  the  correlation  between 
traces  decreases  as  ka  approaches  one.  Also,  because  many  different  wave  types  exist  for 
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these  marine  models,  the  correlation  between  traces  is  range  dependent,  even  for  the 
laterally  homogeneous  case. 
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INTRODUCTION 


Marine  seismic  studies  have  been  paramount  in  unravelling  the  gross  structure  of 
the  oceanic  crust.  Various  models  of  this  structure  have  evolved  and  are  constantly  being 
refined  as  data  collection,  analysis,  and  modeling  techniques  become  more  effective. 
Currently  under  debate  is  the  natrne  of  the  small-scale  lateral  heterogeneity  of  the  crust  as 
well  as  the  depth  dependence  of  the  lateral  heterogeneity.  Li  the  high  frequency  limit,  that 
is,  when  the  structure  is  much  larger  than  the  seismic  wavelength  (  X  ),  seismic  energy 
propagates  coherently  through  the  crust  and  produces  coherent  seismogram  arrivals. 
Information  about  large-scale  features  of  crustal  structure  is  obtained  from  the  arrival  times 
and  slopes  of  these  principal  coherent  arrivals.  However,  as  we  look  for  finer  scale 
structure  it  becomes  apparent  that  this  simple  analysis  of  seismograms  is  no  longer 
adequate. 

Aki  (1982)  defines  scattering  as  "A  process  which  generates  incoherent  signals  by 
three-dimensional  rough  heterogeneities".  As  ka  approaches  unity  (k=2rc/l, 
a=heterogeneity  length  scale)  the  wave  characteristics  of  the  seismic  energy  come  into  the 
picture  causing  interference,  diffractions,  and  in  general  much  more  complex  patterns  of 
propagation.  Crustal  heterogeneities  act  as  scatterers  of  seismic  energy  and  cause  a 
degradation  in  the  coherence  of  propagation  through  the  crust.  Thus,  in  order  to  examine 
the  small  scale  heterogeneity  of  the  crust  it  is  the  incoherence  of  the  expected  arrivals  and 
the  coda  (or  signal  generated  noise)  which  must  be  studied. 

In  general,  most  studies  of  scattering  have  been  carried  out  to  explain  the  amplitude 
attenuation  and  coda  wave  excitation  of  seismic  energy  from  earthquake  or  teleseismic 
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explosive  sources  (Aki,  1973,1982;  Sato,  1982;  McLaughlin,  1983,  McLaughlin 
et.al.,1985,  Frankel  and  Clayton,  1986).  Large  scale  heterogeneities  have  the  greatest 
influence  on  primary  arrival  amplitude  attenuation  while  smaller  scale  heterogeneities  seem 
to  control  coda  wave  frequency  content  and  decay  (Aki,  1973,1982).  Most  analytical  work 
done  to  explain  these  results  depends  on  the  validity  of  the  Bom  approximation  which  only 
considers  the  first  interaction  of  energy  with  a  scatterer  (no  multiple  scattering)  and  that  the 
scattering  is  weak  with  respect  to  the  primary  wave  amplitude  (Hudson  and  Heritage, 
1981). 

These  previous  studies  are  important  foundations  for  this  study  but  are  not  directly 
applicable  for  a  number  of  reasons.  Most  estimates  made  of  heterogeneities  concern  the 
entire  lithosphere  and  parts  of  the  upper  mantle  instead  of  just  the  crust  (upper  5  km  for 
oceanic  crust).  When  the  experiment  is  scaled  to  that  of  typical  marine  refraction  work,  the 
Bom  approximation  is  not  necessarily  valid  (Hudson  and  Heritage,  1981).  The  case  of 
scattering  of  vector  (elastic)  waves  from  a  point  source  travelling  through  a  medium  with 
velocity  gradients,  topography,  and  random  heterogeneities  has  little  chance  to  be  solved 
analytically.  We  wish  to  investigate  heterogeneities  with  length  scales  less  than  the  ranges 
of  typical  seismic  refraction  surveys  (ranges  from  5-50  km). 

An  example  of  the  influence  of  this  size  of  lateral  heterogeneities  on  the  seismic 
signal  is  shown  in  figure  2.1.  The  data  in  figure  2.1  are  from  a  borehole  seismic 
experiment  at  the  Deep  Sea  Drilling  Project  (DSDP)  site  504B  near  the  Costa  Rica  Rift 
(Stephen,  1987).  This  line  of  data  was  collected  using  a  borehole  receiver  placed  within 
the  hole  42  meters  into  basement.  Shots  were  fired  in  a  circle  of  radius  6.0  kilometers  and 
at  azimuths  from  0  to  360  degrees.  Under  the  assumptions  of  lateral  homogeneity  and 
azimuthal  isotropy,  the  arrivals  in  figure  2. 1  should  be  perfecdy  coherent  and  have  exactly 
the  same  arrival  time.  However,  this  is  obviously  not  the  case  for  the  data  shown  here.  In 


Figure  2.1.  Borehole  seismic  data  from  the  Deep  Sea  Drilling  Project  hole  504B  near  the 
Costa  Rica  Rift.  These  data  were  produced  by  placing  a  borehole  receiver  42  meters  into 
basement  and  firing  a  circle  of  surface  shots  (at  azimuths  from  0  to  360  degreees)  at  a 
range  of  6.0  kilometers  from  the  borehole.  The  effects  of  small  scale  lateral  heterogeneity 
in  the  upper  oceanic  crust  appear  as  coda  (signal  generated  noise)  and  azimuth  dependent 
amplitudes,  travel  times,  and  correlation,  (figure  from  Stephen,  1987) 
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fact,  there  is  a  significant  amount  of  signal  generated  noise,  as  well  as  large  amplitude  and 
travel  time  anomalies.  It  is  not  the  intent  of  this  study  to  reproduce  the  data  seen  in  figure 
2.1.  but  rather,  we  are  interested  in  the  types  of  seconday  characteristics  seen  in  these  data. 

In  this  study,  we  will  use  the  finite  difference  method  to  solve  the  forward  problem 
of  seismic  wave  propagation  through  oceanic  crust  with  random  heterogeneity.  Synthetic 
seismograms  calculated  by  this  method  can  be  used  to  quantify  the  change  in  seismic  traces 
due  to  the  heterogeneities.  An  elastic  finite  difference  formulation  of  the  heterogeneous 
wave  equation  (P  and  S  velocity  as  well  as  density  can  be  varied)  allows  us  to  create  any 
realistic  crustal  model  without  the  simplifying  assumptions  of  most  theoretical  scattering 
analyses.  Numerical  Schlieren  diagrams,  or  snapshots,  of  the  energy  propagating  through 
the  models  show  that  a  significant  amount  of  energy  is  backscattered  and  converted  by  the 
heterogeneous  crust.  This  scattered  energy  shows  up  as  'noise'  on  the  synthetic 
seismogram  traces.  The  effect  of  wave  propagation  through  media  with  different  length 
scales  of  heterogeneity  are  presented. 
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FINITE  DIFFERENCE  MODELING 


The  finite  difference  method  provides  an  excellent  solution  to  the  problem  of  wave 
propagation  through  random  media.  The  finite  difference  scheme  used  in  this  paper  is 
based  on  the  scheme  of  Virieux  (1986)  and  includes  compressional  and  shear  velocity 
variations  and  density  variations.  It  was  originally  presented  as  a  solution  to  the  first  order 
system  in  terms  of  particle  velocities  and  stresses,  but  in  order  to  reduce  memory  storage 
requirements  we  rewrote  the  formulation  for  the  second  order  system  in  terms  of  particle 
displacements.  The  calculations  are  carried  out  in  two  dimensions  for  a  line  source  in 
Cartesian  coordinates.  The  equations  solved  are; 
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where  u,w  are  the  horizontal  and  vertical  displacements  respectively,  p  is  density,  x  is 
stress,  and  X,ji  are  Lame's  elastic  parameters.  Centered  finite  differences  are  used  to 
approximate  each  derivative  in  space  and  time.  Dependent  and  independent  variables  are 


101 


m-1 


n-1/2  V 


n+1/2V 


m 


A 


A 


P—Y--* 


A 


A 


*-4-4 


■A- 


A 


m-1/2 


m+1/2 


m+1 

n-1 


V 


n 


V 


n+1 


•  u  ,  p 
■  w,  p 

A  Txx  -  Yz  ’  X+2^’  * 


V  t 


xz 


,  y. 


Figure  2.2.  Location  of  dependent  and  independent  variables  on  the  finite  difference  grid. 
Horizontal  and  vertical  displacements  and  stresses  are  defined  at  different  points  of  the 
grid,  (after  Virieux,  1986) 
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defined  on  the  grid  as  shown  in  figure  2.2.  This  formulation  implicitly  calculates  P  and  SV 
propagation  (SH  is  not  present  in  the  x-z  plane),  compressional  to  shear  conversion, 
interference  effects  and  interface  phenomena. 

We  prefer  this  modified  Virieux  scheme  to  the  Nicoletis  scheme  used  in  previous 
studies  (Nicoletis,  1981,  Bhasavanija,  1983,  Dougherty  and  Stephen,  1987  (chapter  1)) 
because  is  generates  less  parasidc  noise,  particularly  in  the  water  column.  This  noise,  for 
the  Nicoletis  scheme  is  not  large  enough  to  contaminate  time  series  results  but  it  does  spoil 
the  appearance  of  snapshots  and  creates  non-zero  curl  of  the  displacement  field  in  the  water 
column.  These  issues  are  discussed  further  in  Stephen  (1988).  A  complete  comparison  of 
this  scheme  with  fourth  order  schemes  (Alford  et.al.,  1974,  Frankel  and  Clayton,  1986) 
has  not  been  made  by  these  authors  but  Fomberg  (1987)  has  given  some  results  which 
suggest  that  the  fourth-order  scheme  has  advantages.  Since  Fomberg  (1987)  did  not 
investigate  staggered  grids,  as  used  here,  his  results  should  be  considered  preliminary. 
Computer  time  and  space  requirements  for  the  Virieux  scheme  were  not  prohibitive  for 
these  models  and  the  accuracy  of  the  seismograms  for  the  Virieux  method  are  excellent.  A 
complete  discusion  of  these  issues  is  beyond  the  scope  of  this  paper. 

Grid  layout  and  boundary  conditions  are  given  in  figure  2.3.  The  grid  is  split  into 
two  homogeneous  layers  and  a  heterogeneous  transition  zone.  Constant  parameters  in  the 
homogeneous  water  (Vp=1.5  km/s,  Vs=0.0,  p=1.0)  and  deep  crustal  (Vp=6.0,  Vs=3.46, 
p=2.525)  layers  allow  a  simplification  of  equation  1  and  a  saving  in  computation  time. 
Compressional  and  shear  velocities  as  well  as  density  can  be  varied  at  each  grid  point 
within  the  transition  zone. 

The  initial  conditions  are  zero  particle  displacement  and  velocity  everywhere  on  the 
grid.  Boundary  conditions  for  the  models  were  chosen  to  minimize  numerical  aritifacts 
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returning  from  the  edges  of  the  grid.  The  right  hand  side  of  the  grid  is  an  absorbing 
boundary  based  on  paraxial  approximations  of  the  wave  equation  given  by  Clayton  and 
Engquist  (1977).  This  type  of  boundary  was  found  to  be  unstable  for  the  bottom 
boundary  of  the  grid  at  large  times.  Therefore,  reflections  from  the  bottom  boundary  are 
suppressed  by  using  a  form  of  the  telegraph  equation  combined  with  the  homogeneous 
elastic  wave  equation  (Levander,  1985,  Ceijan  et.al.,  1985).  Two  attenuation  terms  are 
added  to  the  wave  equation  which  successively  decrease  wave  amplitude  in  a  region  near 
the  bottom  boundary.  The  section  of  the  grid  which  includes  the  telegraph  formulation  is 
shown  by  the  stipled  area  in  figure  2.3. 

The  effect  of  the  line  source  is  introduced  as  a  time  dependent  boundary  condition 
on  the  top  of  the  grid.  A  second  grid  is  introduced  along  the  top  edge  (see  figure  2.3)  to 
absorb  upcoming  energy  and  to  prevent  it  from  reflecting  back  into  the  region  of  interest  as 
a  false  arrival.  Up  and  downgoing  wavefields  are  separated  as  in  Alterman  and  Lowenthal 
(1972) .  The  upgoing  field  is  then  absorbed  at  the  top  of  the  second  grid  using  the  Clayton 
and  Engquist  (1977)  absorbing  boundary  condition. 

As  stated  above,  a  source  function  is  introduced  at  the  top  of  the  grid,  along  a  layer 
0.5  km  (50  grid  points  at  10  meters  per  grid  point)  above  the  seafloor  (see  figure  2.3).  The 
source  waveform  simulates  a  10  Hz  shot  at  the  sea  surface  and  is  shown  in  figure  2.4.  The 
shape  of  the  pressure  source  wave  is  based  on  the  second  derivative  of  a  Gaussian  shape; 

(2)  p  (x,z,t)  =  ( -A  /  (  4rta2R  ) )  g"(  t  -  R/a  ) 


where 
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g”(  t )  =  4^2  (  3T  -  2$T3  )  exp  ( -$T2  ) 


R  =  (  x5  +  z2  )  J/2  ,  T  =  t  -  ts 

In  equation  2,  p  is  the  pressure  field  at  time  t  and  coordinates  x,z,  A  is  a  unit  constant  with 
dimensions  (mass  x  length2  /  time),  a  is  the  compressional  velocity  of  water,  g  is  a 
Gaussian  curve,  ts  is  a  time  shift  parameter,  and  £  is  a  pulse  width  parameter  for  frequency 
scaling  (Stephen  et.al.,  1985,  Kelly  et.al.,  1979).  A  value  of  657  for  £  was  used  which 
produces  a  signal  with  peak  frequency  of  10  Hz  and  an  upper  half  power  frequency  of  1 3.5 
Hz.  The  time  increment  used  was  0.001  seconds. 

A  number  of  different  output  formats  are  available  from  the  finite  difference 
method.  In  fact,  one  of  the  advantages  of  the  method  is  the  ability  to  save  panicle 
displacements  for  any  number  of  points  on  the  grid  and  for  all  times.  A  particularly 
informative  way  to  view  wave  propagation  through  the  models  is  by  using  numerical 
Schlieren  diagrams.  These  'snapshots'  of  the  wave  field  travelling  through  the  grid  are 
produced  by  saving  the  value  of  horizontal  and  vertical  displacement  at  each  point  in  the 
grid  for  a  given  time  step.  The  quantities  plotted  in  the  Schlieren  diagrams  are  related  to 
compressional  and  shear  energies  and  are  calculated  by  using  the  spatial  divergence  and 
curl  operators  on  the  displacements.  Morse  and  Feshbach  (1953)  define  compressional 
(Ec)  and  shear  (Es)  energy  as; 

(4)  Ec  =  (  X+2 4  )  (  V.u  )  2 

(5)  ES  =  M-Vxu)2 
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Figure  2.3.  Layout  of  the  finite  difference  grid  for  the  models  presented.  The  source  is 
introduced  along  a  layer  50  grid  points  above  the  water-solid  interface.  Homogeneous 
layers  of  water  (top  of  grid)  and  crust  (bottom)  are  solved  using  a  formulation  of  the 
homogeneous  wave  equation.  The  heterogeneous  transition  zone  contains  all  lateral  and 
vertical  velocity  (both  compressional  and  shear)  and  density  variations.  Energy  is 
attenuated  by  using  the  telegraph  equation  in  the  stipled  area  at  the  bottom  of  the  grid. 
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Figure  2.4.  Pressure  source  function  used  in  the  finite  difference  models. 
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In  order  to  preserve  the  divergence  and  curl  sign  information  while  showing  relative 
compressional  and  shear  energy  amplitudes,  we  plot  the  following  normalized  quantities  in 
the  numerical  Schlieren  diagrams; 


(6) 


Norm  P  =  sign  (  V»u  )  *  (  Ec ) 


1/2 


(7) 


Norm  S  =  sign  (  Vxu  )  *  (  Es  ) 


1/2 


Synthetic  seismograms  are  produced  by  saving  displacements  and/or  pressure  at 
any  point  on  the  grid  for  a  series  of  timesteps.  Time  series  of  horizontal  and  vertical 
displacements,  as  well  as  pressure  were  produced  for  a  row  of  'receivers'  20  meters  above 
the  water-solid  interface  and  a  row  30  meters  below  the  interface.  This  was  done  to 
compare  the  pressure  signal  at  a  water  column  hydrophone  to  displacements  from  a  buried 
geophone.  Receivers  were  placed  every  four  grid  points  (equivalent  to  40  meter  spacing) 
for  each  row  and  values  were  saved  for  every  fourth  timestep  (.004  second  sampling 
interval). 

SCATTERING  MODELS 


A  number  of  models  were  run  with  different  length  scales  of  heterogeneity  and 
spatial  filter  types  to  examine  the  effect  on  the  seismograms.  A  control  model,  FLAT,  was 
first  run  to  demonstrate  wave  propagation  through  an  upper  oceanic  crust  with  no  lateral 
heterogeneity.  Figure  2.5  presents  the  velocity-depth  function  used  for  FLAT.  The 
velocity-depth  function  is  based  on  results  of  Purdy  (1982)  and  Ewing  and  Purdy  (1982) 
for  young  ocean  crust  studied  in  the  Rivera  Ocean  Seismic  Experiment.  Dougherty  and 
Stephen  (1987,  chapter  1)  used  a  similar  V-Z  function  to  model  deterministic  scattering 
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from  seafloor  features.  All  of  the  models  assume  a  Poisson's  ratio  of  0.25  and  the  density 
(  r  )  is  calculated  from  P-wave  velocity  by  the  following  relation  (Nafe  and  Drake,  1957); 


(5) 


p  (  x,z  )  =  . 252  +  (. 3788  )*  Vp(x,z) 


Lateral  heterogeneity  was  created  by  adding  a  series  of  normally  distributed  random 
velocity  variations  to  the  V-Z  function  of  FLAT.  Each  grid  point  in  the  crustal  section  of 
the  heterogeneous  transition  zone  has  a  random  component  added  to  it  The  random 
velocity  variations  have  zero  mean  so  that  the  average  velocity  for  any  given  layer  is 
approximately  equal  to  the  velocity  of  that  layer  in  the  laterally  homogeneous  case.  The 
heterogeneity  standard  deviation  is  +/-  10  percent  at  the  water  solid  interface  and  decreases 
down  to  zero  at  a  depth  of  1.3  km.  (130  grid  points )  below  the  interface.  The  decrease  in 
heterogeneity  with  depth  is  meant  to  simulate  the  sealing  of  cracks  and  fissures  with  depth 
in  the  upper  oceanic  crust. 

After  adding  random  velocity  variations  to  each  grid  point  in  the  transition  zone,  the 
entire  velocity  field  was  spatially  filtered  to  create  different  heterogeneity  correlation 
lengths.  A  filter  which  produced  a  Gaussian  autocorrelation  function  was  used  for  most  of 
the  models  in  this  study.  Spatial  correlation  lengths  of  50, 100,  and  200  meters  where 
used  for  models  GAUSS50,  GAUSS  100,  and  GAUSS200  respectively.  One  model  with 
a  self-similar  autocorrelation  function  and  a  correlation  length  of  200  meters  is  also 
presented  (SELFSIM200).  The  self-similar  distribution  is  based  on  the  fractal  geometry  of 
Mandelbrot  (1977)  and  contains  uniform  fluctuations  over  a  range  of  length  scales. 

Frankel  and  Clayton  (1986)  describe  these  filter  correlation  functions  in  terms  of  their 
definitions  and  fluctuation  power  spectra. 
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Figure  2.6  shows  the  velocity  fields  of  the  five  heterogeneous  models.  GAUSS  10 
has  just  the  original  unfiltered  random  number  sequence  added  to  the  gradients  of  FLAT. 
Although  this  would  be  an  uncorrelated  sequence  for  the  continuous  case  ,  the  discrete  1 0 
meter  grid  spacing  of  the  finite  difference  formulation  imposes  a  correlation  length  of  10 
meters  on  the  velocity  field  of  GAUSS  10.  GAUSS50,  GAUSS100,  and  GAUSS200  are 
based  on  the  same  initial  random  number  sequence  but  are  filtered  with  different  correlation 
lengths  of  a  Gaussian  filter.  SELFSIM200  has  a  self-similar  aoutocorrelation  function  with 
a  correlation  length  of  200  meters. 

Laterally  homogeneous  model 

Before  examining  propagation  through  the  random  models,  it  is  useful  to  review 
energy  propagation  and  partitioning  for  this  seafloor  geometry  with  the  laterally 
homogeneous  case  (FLAT).  Unlike  propagation  of  plane  or  spherical  waves  through  a 
halfspace,  simulations  of  marine  refraction  experiments  contain  a  water-solid  interface 
underlain  by  a  region  of  velocity  gradient.  This  drastic  change  in  elastic  parameters  results 
in  a  number  of  interesting  partitioning,  conversion,  and  interface  effects.  Figure  2.7  is  a 
series  of  numerical  Schlieren  diagrams  which  show  the  propagation  of  normalized 
cr  mpressional  and  shear  energy  across  the  FLAT  velocity  field.  These  'snapshots'  are 
shown  for  0.4  second  intervals  starting  at  2.0  seconds  after  the  shot.  The  top  of  each 
frame  corresponds  to  a  depth  of  2.5  km  below  the  sea  surface  with  the  sea  floor  at  3.0  km. 
depth.  The  bottom  of  the  grid  corresponds  to  a  total  depth  of  5.5  km.  (2.5  km  below  the 
water-solid  interface). 


Ill 


Figure  2.5.  Compressional  velocity  field  for  the  laterally  homogeneous  model  FLAT.  P- 
velocity  (a)  at  the  water-solid  interface  is  3.2  km/sec  and  increases  to  6.0  km/sec  at  a  depth 
of  1.3  km  below  the  interface  (total  depth  4.3  km).  A  steep  gradient  of  4.0  sec*l  is  present 
from  the  interface  to  0.4  km  below  the  interface,  with  a  second  less  steep  gradient  of  1.33 
sec'l  from  0.4  to  1.3  km.  below  the  interface.  It  is  assumed  that  the  crust  is  a  Poisson's 
solid  so  that  the  shear  velocity,  ^=0/(3)^, 
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Figure  2.6.  Compressional  velocity  fields  for  the  random  heterogeneity  models 
GAUSS  10-GAUSS200.  In  each  case,  the  average  velocity  for  each  depth  is  equal  to  the 
velocity  at  that  depth  in  the  laterally  homogeneous  model  FLAT  (figure  2.5).  Also,  the 
standard  deviation  of  the  random  velocities  is  +/-  10  percent  at  the  water-solid  interface, 
decreasing  to  zero  at  a  depth  of  4.3  km  (1.3  km  below  the  interface).  a.GAUSSlO, 
Gaussian  random  perturbations  applied  to  each  grid  point,  since  no  spatial  filtering  is  done, 
the  effective  correlation  length  for  GAUSS  10  is  10  meters  (equal  to  the  grid  spacing),  b. 
GAUSS50,  velocity  field  of  GAUSS  10  filtered  with  a  2-D  Gaussian  spatial  filter, 
correlation  length  equal  to  50  meters,  c.  GAUSS  100,  Gaussian  spatial  filter,  correlation 
length  of  100  meters,  d.  GAUSS200,  Gaussian  spatial  filter,  correlation  length  200 
meters,  e.  SELFSIM200,  self-similar  spatial  filter,  correlation  length  200  meters.  Note  in 
the  latter  case,  the  presence  of  smaller  scale  roughness  superimposed  on  the  larger  features 
of  the  filtered  field. 
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Figure  2.6.  Compressional  velocity  fields  for  the  random  heterogeneity  models 
GAUSS  10-GAUSS200.  In  each  case,  the  average  velocity  for  each  depth  is  equal  to  the 
velocity  at  that  depth  in  the  laterally  homogeneous  model  FLAT  (figure  2.5).  Also,  the 
standard  deviation  of  the  random  velocities  is  +/-  10  percent  at  the  water-solid  interface, 
decreasing  to  zero  at  a  depth  of  4.3  km  (1.3  km  below  the  interface).  a.GAUSSlO, 
Gaussian  random  perturbations  applied  to  each  grid  point,  since  no  spatial  filtering  is  done, 
the  effective  correlation  length  for  GAUSS  10  is  10  meters  (equal  to  the  grid  spacing),  b. 
GAUSS50,  velocity  field  of  GAUSS  10  filtered  with  a  2-D  Gaussian  spatial  filter, 
correlation  length  equal  to  50  meters,  c.  GAUSS  100,  Gaussian  spatial  filter,  correlation 
length  of  100  meters,  d.  GAUSS200,  Gaussian  spatial  filter,  correlation  length  200 
meters,  e.  SELFSIM200,  self-similar  spatial  filter,  correlation  length  200  meters.  Note  in 
the  latter  case,  the  presence  of  smaller  scale  roughness  superimposed  on  the  larger  features 
of  the  filleicd  Held. 
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FLAT  represents  a  section  of 'typical'  young  ocean  crust  with  no  lateral 
heterogeneity.  With  a  seafloor  at  3  km  depth,  it  takes  2  seconds  for  the  source  to  reach  the 
seafloor,  as  can  be  seen  in  figure  2.7a.  After  incidence  on  the  water-solid  interface,  the 
direct  compressional  wave  from  the  source  splits  into  a  direct  wave,  a  seafloor  reflected  P- 
wave,  a  transmitted  P-wave,  and  a  converted  transmitted  shear  wave  (figure  2.7b).  The 
transmitted  waves  'dive'  in  the  presence  of  velocity  gradients  and  are  eventually  turned 
upwards  by  refraction  to  be  incident  on  the  seafoor  from  below.  These  waves  are  referred 
to  as  P  and  S  diving  waves.  After  2.8  seconds  (figure  2.7c)  some  of  the  P-diving  energy 
has  returned  to  the  seafloor  from  below  causing  two  effects.  First,  a  plane  wave  is  set  up 
in  the  water  from  compressional  energy  transmission  through  the  crust.  Second  there  is  a 
P-S  converted  wave  caused  by  the  P-diving  wave  interaction  with  the  seafloor.  This  wave 
is  the  P1P2S2  wave  (Brekhovskikh,  1960)  and  is  seen  coupled  to  the  P-diving  wave  at  the 
water-solid  interface  and  tangential  to  the  S-diving  wave  at  depth.  Of  course,  there  is  also  a 
downward  P-reflection  from  the  P-diving  interaction  with  the  seafloor,  but  at  these  early 
times  the  P-diving  multiples  are  indistinguishable  from  the  initial  P-arrival  (Cerveny  and 
Ravindra,  1971,  Dougherty  and  Stephen,  1987).  The  emergence  of  the  P-interference  head 
wave  which  results  from  these  multiples  can  be  seen  in  the  3.2  second  snapshot  (figure 
2.7d). 


Eventually,  the  S-diving  wave  is  turned  back  to  be  incident  on  the  seafloor  from 
below  and  shows  the  same  general  characteristics  of  the  P-diving  wave  (see  figures  2.7d- 
h).  At  post  critical  ranges  (that  is,  after  the  P  and  S-diving  waves  have  separated  from  the 
direct  wave),  an  evanescent  direct  wave  root  is  established  which  is  coupled  to  the  direct 
wave  at  the  seafloor  (figures  2.7e-h).  This  wave  has  both  compressional  and  shear 
components  and  dies  off  exponentially  with  depth  (Stephen  and  Bolmer,  1985).  The  only 
other  major  wave  to  develop  is  the  pseudo- Rayleigh  wave  which  appears  in  the  4.4  and  4.8 
second  snapshots  (figures  2.7g-h). 
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Synthetic  seismograms  from  model  FLAT  are  shown  in  figure  2.8  for  hydrophones 


and  buried  receivers  along  the  seafloor.  The  Schlieren  diagrams  for  FLAT  reveal  that  the 
interactions  between  the  different  waves  present  is  quite  complex  even  for  this  laterally 
homogeneous  case.  However,  it  is  important  to  note  that  all  of  the  waves  appear  as 
coherent  arrivals  in  the  seismograms.  Deterministic  scattering  from  laterally  homogeneous 
crust  is  complex  but  predictable.  This  is  evidenced  by  apparently  strong  correlation 
between  seismic  traces  in  figure  2.8. 


Random  models 

The  object  of  this  study  is  to  demonstrate  the  degradation  in  the  seismic  signal 
which  is  caused  by  wave  scattering  from  random  velocity  variations  in  the  upper  oceanic 
crust.  A  series  of  models  with  Gaussian  velocity  autocorrelation  functions  were  ran  to 
demonstrate  the  effect  on  the  seismic  signal.  The  Gaussian  autocorrelation  function  was 
chosen  for  most  of  the  models  because  it  is  a  particularly  well  studied  problem  (Chernov, 
1960).  One  model  with  a  self-similar  (fractal)  distribution  of  velocity  fluctuations  is  also 
presented. 

Model  GAUSS50,  as  with  all  of  the  following  models,  contains  a  +/-  10  percent 
random  velocity  fluctuation  at  the  water-solid  interface.  The  fluctuation  decreases  to  zero  at 
a  depth  of  1.3  km.  below  the  seafloor  (see  figure  2.6).  The  velocity  field  of  GAUSS50 
has  a  Gaussian  autocorrelation  function  with  a  characteristic  correlation  length  of  50 
meters.  Since  a  predominantly  10  Hz  source  is  used  in  these  models,  ka  is  approximately 
equal  to  one  at  the  seafloor.  Wavefront  snapshots  for  model  GAUSS50  are  presented  in 
figure  2.9.  Synthetic  seismograms  for  water  column  and  buried  receivers  are  shown  in 
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Figure  2.7.  Numerical  Schlieren  diagrams  of  normalized  compressional  and  shear  (top 
and  bottom,  respectively)  energy  for  the  laterally  homogeneous  model  FLAT.  The  frames 
correspond  to  the  ranges  0-10  km.  in  the  model  geometry  shown  in  figure  2.3.  After  a 
simulated  surface  shot,  the  direct  wave  is  just  reaching  the  seafloor  after  2.0  seconds  (fig. 
2.7a)  and  partitions  into  direct,  reflected,  transmitted  P  and  transmitted  S  waves  after 
interaction  with  the  seafloor  (fig.  2.7b).  The  P1P2S2  arrival  is  seen  after  2.8  seconds  (fig. 
2.7c)  and  the  interference  P-wave  is  starting  to  be  distinct  after  3.2  seconds  (fig.  2.7d).  P 
and  S  components  of  the  direct  wave  root  show  up  as  evanescent  waves  coupled  to  the 
direct  wave  after  the  P  and  S  diving  waves  have  separated  from  the  direct  wave  (fig.  2.7f- 
h).  The  only  other  major  wave  to  occur  is  the  pseudo-Rayleigh  wave  seen  in  figure  2.7g- 
h. 
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Figure  2.7.  Numerical  Schlieren  diagrams  of  normalized  compressional  and  shear  (top 
and  bottom,  respectively)  energy  for  the  laterally  homogeneous  model  FLAT.  The  frames 
correspond  to  the  ranges  0-10  km.  in  the  model  geometry  shown  in  Figure  2.3.  After  a 
simulated  surface  shot,  the  direct  wave  is  just  reaching  the  seafloor  after  2.0  seconds  (Fig. 
2.7a)  and  partitions  into  direct,  reflected,  transmitted  P  and  transmitted  S  waves  after 
interaction  with  the  seafloor  (Fig.  2.7b).  The  P1P2S2  arrival  is  seen  after  2.8  seconds  (Fig. 
2.7c)  and  the  interference  P-wave  is  starting  to  be  distinct  after  3.2  seconds  (Fig.  2.7d).  P 
and  S  components  of  the  direct  wave  root  show  up  as  evanescent  waves  coupled  to  the 
direct  wave  after  the  P  and  S  diving  waves  have  separated  from  the  direct  wave  (Fig.  2.7f- 
h).  The  only  other  major  wave  to  occur  is  the  pseudo-Rayleigh  wave  seen  in  Figure  2.7g- 
h. 
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Figure  2.7.  Numerical  Schlieren  diagrams  of  normalized  compressional  and  shear  (top 
and  bottom,  respectively)  energy  for  the  laterally  homogeneous  model  FLAT.  The  frames 
correspond  to  the  ranges  0-10  km.  in  the  model  geometry  shown  in  figure  2.3.  After  a 
simulated  surface  shot,  the  direct  wave  is  just  reaching  the  seafloor  after  2.0  seconds  (fig. 
7a)  and  partitions  into  direct,  reflected,  transmitted  P  and  transmitted  S  waves  after 
interaction  with  the  seafloor  (fig.  7b).  The  P1P2S2  arrival  is  seen  after  2.8  seconds  (fig. 
7c)  and  the  interference  P-wave  is  starting  to  be  distinct  after  3.2  seconds  (fig.  7d).  P  and 
S  components  of  the  direct  wave  root  show  up  as  evanescent  waves  coupled  to  the  direct 
wave  after  the  P  and  S  diving  waves  have  separated  from  the  direct  wave  (fig.  7f-h).  The 
only  other  major  wave  to  occur  is  the  pseudo-Rayleigh  wave  seen  in  figure  7g-h. 
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Figure  2.7.  Numerical  Schlieren  diagrams  of  normalized  compressional  and  shear  (top 
and  bottom,  respectively)  energy  for  the  laterally  homogeneous  model  FLAT.  The  frames 
correspond  to  the  ranges  0-10  km.  in  the  model  geometry  shown  in  figure  2.3.  After  a 
simulated  surface  shot,  the  direct  wave  is  just  reaching  the  seafloor  after  2.0  seconds  (fig. 
2.7a)  and  partitions  into  direct,  reflected,  transmitted  P  and  transmitted  S  waves  after 
interaction  with  the  seafloor  (fig.  2.7b).  The  P1P2S2  arrival  is  seen  after  2.8  seconds  (fig. 
2.7c)  and  the  interference  P-wave  is  starting  to  be  distinct  after  3.2  seconds  (fig.  2.7d).  P 
and  S  components  of  the  direct  wave  root  show  up  as  evanescent  waves  coupled  to  the 
direct  wave  after  the  P  and  S  diving  waves  have  separated  from  the  direct  wave  (fig.  2.7f- 
h).  The  only  other  major  wave  to  occur  is  the  pseudo-Rayleigh  wave  seen  in  figure  2.7g- 
h. 
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Figure  2.8.  Synthetic  seismograms  for  model  FLAT.  a.  Pressure  signal  for  receivers  20 
meters  above  the  seafloor,  b.  Vertical  displacement  for  receivers  30  meters  below  the 
seafloor,  c.  Horizontal  displacement  for  receivers  30  meters  below  the  seafloor. 
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figure  2.10. 


An  examination  of  figures  2.9  and  2.10  reveals  that  the  principal  arrivals  present  in 
the  laterally  homogenous  case  still  dominate  the  results.  The  amplitude  and  coherence  of 
the  P  and  S  diving  waves  have  decreased  due  to  the  random  heterogeneities.  However,  the 
amount  of  lateral  heterogeneity  is  not  great  enough  to  significantly  alter  the  arrival  times  or 
slopes  of  the  major  arrivals  in  the  seismograms.  One  arrival  which  does  seem  to  be  greatly 
affected  is  the  interference  P-head  wave.  The  laterally  heterogeneous  velocity  field  of 
GAUSS50  causes  a  more  destructive  interference  between  the  diving  wave  multiples  and 
the  interference  head  waves  are  not  as  strong  as  in  the  laterally  homogeneous  case. 
Dougherty  and  Stephen  (1987,  chapter  1)  also  observed  this  destructive  interference  for  the 
case  of  simple  seafloor  topographic  scatterers.  This  explains  why  interference  head  waves 
have  not  been  identified  in  field  data. 

Random  velocity  fluctuations  of  the  magnitude  and  spatial  dimensions  of 
GAUSS50  scatter  seismic  energy  into  body  waves  (both  compressional  and  shear  types) 
and  interface  waves.  Body  wave  scattering  of  compressional  and  shear  energy  is  the 
dominant  secondary  feature  in  the  snapshots  of  figure  2.9  (after  the  expected  primary 
arrivals).  A  complex  pattern  of  backscattered,  multiply  scattered,  and  converted  energy  is 
produced  between  the  principle  wavefronts.  Converted  shear  energy  is  found  before  the  S- 
transmitted  wave  from  scattering  of  the  P-diving  wave  into  both  P  and  S  scattered  waves. 
Although  much  of  the  body  wave  scattering  is  incoherent,  some  faint  coherent  waves  can 
be  seen  cylindrically  spreading  below  the  seafloor  in  the  later  timesteps  (figures  2.9c-h).  It 
appears  that  most  of  the  energy  scattered  into  body  waves  is  shear  energy,  either  from  a 
scattered  shear  arrival  or  from  conversion  during  P-wave  scattering. 
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A  great  deal  of  energy  is  scattered  by  the  lateral  velocity  fluctuations  back  into  the 
water  column.  The  backscattered  body  waves  appear  as  plane  waves  in  the  water  column 
and  show  up  in  the  time  series  as  relatively  incoherent  noise’.  Upon  closer  examination  of 
the  time  series,  coherent  sections  of  scattered  body  waves  can  be  observed.  Most  of  the 
coherent  scattered  body  waves  are  shear  arrivals  and  travel  with  a  speed  of  approximately 
2.8  km/sec.  Figure  2.1 1  is  a  plot  of  two  sections  of  the  time  series  from  receivers  with  a 
spacing  of  40  meters  (instead  of  480  meter  spacing  in  Figures  2.8  and  2.10).  Small 
sections  of  coherent  S-phases  are  present  which  are  not  obviously  coherent  when  observed 
with  a  larger  seismogram  spacing. 

Another  effect  of  secondary  scattering  of  the  primary  arrivals  is  the  generation  of 
secondary  Stoneley  waves  along  the  water-solid  interface.  These  waves  do  not  appear  in 
the  laterally  homogeneous  case  because  the  only  source  is  located  too  far  from  the  interface. 
However,  in  the  laterally  heterogeneous  models,  velocity  fluctuations  of  a  significant 
magnitude  and  spatial  dimension  can  act  as  secondary  sources  of  seismic  energy  when 
encountered  by  primary  source  energy.  If  these  secondary  sources  are  in  close  enough 
proximity  to  the  water-solid  interface,  secondary  Stoneley  waves  can  be  generated. 

Secondary  Stoneley  waves  in  GAUSS50  appear  to  be  generated  by  diffractions  of 
the  P-diving  wave,  the  direct  wave  root,  and  the  S-diving  wave.  The  amplitude 
characteristics,  velocity,  and  particle  motion  of  Stonely  waves  distinguish  them  from  the 
other  arrivals  present  The  amplitude  of  these  interface  waves  dies  off  exponentially  with 
distance  away  from  both  sides  of  the  water-solid  interface.  Stoneley  waves  appear  as  small 
packets  of  energy  propagating  along  the  interface  rather  than  as  wave  fronts  travelling 
through  the  bodies  of  water  or  solid  in  the  models  (see  figure  2.9). 
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Figure  2.9.  Numerical  Schlieren  diagrams  of  normalized  compressional  and  shear  energy 
for  model  GAUSS50.  The  velocity  field  has  Gaussian  autocorrelated  random  velocity 
variations  with  correlation  length  of  50  meters  (see  figure  2.6  for  velocity  field).  The 
principle  waves  are  very  similar  to  those  seen  in  the  laterally  homogeneous  case  of  figure 
2.7.  Secondary  effects  include  incoherent  backscattering  of  P  and  S  energy  from  the  P  and 
S  diving  waves  and  Stoneley  wave  generation  along  the  interface  from  the  transmitted 
waves  and  the  direct  wave  root.  Note  the  predominance  of  randomly  scattered  shear 
energy  from  both  the  S-diving  wave  and  converted  from  the  P-diving  wave.  Shear  energy 
from  the  P-diving  wave  appears  in  front  of  the  S-diving  wave. 
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Figure  2.9.  Numerical  Schlieren  diagrams  of  normalized  compressional  and  shear  energy 
for  model  GAUSS50.  The  velocity  field  has  Gaussian  autocorTelated  random  velocity 
variations  with  correlation  length  of  50  meters  (see  figure  2.6  for  velocity  field).  The 
principle  waves  are  very  similar  to  those  seen  in  the  laterally  homogeneous  case  of  figure 
2.7.  Secondary  effects  include  incoherent  backscattering  of  P  and  S  energy  from  the  P  and 
S  diving  waves  and  Stoneley  wave  generation  along  the  interface  from  the  transmitted 
waves  and  the  direct  wave  root.  Note  the  predominance  of  randomly  scattered  shear 
energy  from  both  the  S-diving  wave  and  converted  from  the  P-diving  wave.  Shear  energy 
from  the  P-diving  wave  appears  in  front  of  the  S-diving  wave. 
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Figure  2.9.  Numerical  Schlieren  diagrams  of  normalized  compressional  and  shear  energy 
for  model  GAUSS50.  The  velocity  field  has  Gaussian  autocorrelated  random  velocity 
variations  with  correlation  length  of  50  meters  (see  figure  2.6  for  velocity  field).  The 
principle  waves  are  very  similar  to  those  seen  in  the  laterally  homogeneous  case  of  figure 
2.7.  Secondary  effects  include  incoherent  backscattering  of  P  and  S  energy  from  the  P  and 
S  diving  waves  and  Stoneley  wave  generation  along  the  interface  from  the  transmitted 
waves  and  the  direct  wave  root.  Note  the  predominance  of  randomly  scattered  shear 
energy  from  both  the  S -diving  wave  and  converted  from  the  P-diving  wave.  Shear  energy 
from  the  P-diving  wave  appears  in  front  of  the  S-diving  wave. 
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Figure  2.9.  Numerical  Schlieren  diagrams  of  normalized  compressional  and  shear  energy 
for  model  GAUSS50.  The  velocity  field  has  Gaussian  autocorrelated  random  velocity 
variations  with  correlation  length  of  50  meters  (see  figure  2.6  for  velocity  field).  The 
principle  waves  are  very  similar  to  those  seen  in  the  laterally  homogeneous  case  of  figure 
2.7.  Secondary  effects  include  incoherent  backscattering  of  P  and  S  energy  from  the  P  and 
S  diving  waves  and  Stoneley  wave  generation  along  the  interface  from  the  transmitted 
waves  and  the  direct  wave  root.  Note  the  predominance  of  randomly  scattered  shear 
energy  from  both  the  S-diving  wave  and  converted  from  the  P-diving  wave.  Shear  energy 
from  the  P-diving  wave  appears  in  front  of  the  S-diving  wave. 
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The  time  series  in  figure  2. 10  contain  a  number  of  secondary  Stoneley  arrivals. 
Their  low  velocity  and  the  number  present  (especially  after  the  direct  wave)  make  them  hard 
to  identify  with  the  receiver  spacing  shown  in  figure  2.10.  Stoneley  wave  arrivals  become 
more  apparent  when  signals  from  all  of  the  receivers  are  plotted  (40  meter  spacing  as  in 
figure  2.1 1).  The  distinctive  features  of  these  arrivals  are  their  low  velocity  (1.2-1.33 
km/sec), their  linear  moveout  and  their  exponential  decay  in  amplitude  away  from  both  sides 
of  the  interface.  Distinct  Stoneley  wave  arrivals  are  seen  following  the  P-diving  and  direct 
waves  in  blown-up'  sections  of  the  seismograms  (figures  2.1  la-b).  Stoneley  waves 
generated  by  the  P-diving  wave  are  weak  but  significant  when  compared  to  the  P-diving 
amplitude  (figure  2. 1  lb).  A  number  of  Stoneley  waves  are  generated  by  the  direct  wave 
root  and  shear  transmitted  waves  (figure  2. 1  la).  These  are  very  strong  in  amplitude  when 
compared  with  the  P-diving  wave  and  the  Stoneley  waves  generated  by  it.  Although  most 
of  the  Stoneley  waves  are  travelling  to  the  right  across  the  grid,  there  are  a  few  which  are 
scattered  backwards  by  the  direct  wave  root. 

Secondary  Stoneley  waves  have  a  characteristic  elliptical,  generally  retrograde, 
particle  motion  at  the  water-solid  interface  (Bullen  and  Bolt,  1985,  Schirmer,  1980,  Tuthill 
et.al.,  1981).  Elliptical  particle  motion  entails  a  90  degree  phase  shift  between  vertical  and 
horizontal  components.  Thus,  these  arrivals,  although  not  easily  seen  as  coherent  arrivals 
when  the  receiver  spacing  is  too  great,  can  be  identified  by  observing  this  phase  shift  on  a 
superposition  of  vertical  and  horizontal  particle  motion  traces.  Also,  when  plotted  against 
each  other,  the  vertical  and  horizontal  particle  motions  yield  the  charateristic  retrograde 
elliptical  particle  motion  (figure  2. 12). 

The  4.0  second  snapshot  for  the  other  models  with  Gaussian  velocity 
autocorrelation  functions  are  shown  in  figure  2. 13.  The  corresponding  pressure  time  series 
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are  presented  in  figure  2.14.  Again,  as  with  GAUSS50,  all  of  the  expected,  deterministic 
arrivals  occur  as  well  as  a  number  of  secondary  scattered  arrivals.  In  general,  it  appears 
that  the  amount  of  scattered  energy  decreases  as  ka  differs  from  one.  Scattered  arrivals 
become  more  coherent  and  deterministic  as  ka  increases  (ka  appr.  equal  to  4  for 
compressional  waves  in  GAUSS200).  Amplitude  fluctuations  and  residual  travel  time  of 
the  initial  P-diving  wave  also  increase  with  increasing  ka. 

Model  SELFSIM200,  with  a  self-similar  (fractal)  velocity  autocorrelation  function 
is  presented  in  figure  2. 15.  SELFSIM200,  like  GAUSS200  has  a  200  meter  correlation 
length  and  the  large  scale  features  of  the  snapshots  and  time  series  appear  very  similar 
between  the  two.  However,  the  self-similar  distribution  also  contains  smaller  scale 
fluctuations  (see  figure  2.6  for  velocity  fields)  which  scatter  energy  as  in  the  Gaussian 
models  with  smaller  correlation  lengths.  Although  the  principle  effects  of  scattering  seem 
to  be  characteristic  of  the  large  correlation  length  (200m),  some  small  correlation  length 
effects  are  also  present  with  the  self-similar  distribution. 
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Figure  2.10.  Synthetic  seismograms  for  model  GAUSS50.  The  major  arrivals,  as  in  the 
wavefront  snapshots  (fig.  2.9)  are  very  similar  to  those  of  the  laterally  homogeneous  case 
(fig.  2.8).  a.  Pressure  signal  from  receivers  20  meters  above  the  interface,  b.  Vertical 
displacement  from  receivers  30  meters  below  the  interface,  c.  Horizontal  displacement 
from  receivers  30  meters  below  the  interface.  Only  a  few  Stoneley  waves  are  coherent  with 
this  receiver  spacing  (480  meters). 
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Figure  2.1 1.  Closeup  sections  of  the  pressure  synthetic  seismograms  for  model 
GAUSS50.  Receiver  spacing  for  the  seismograms  in  this  figure  is  40  meters  (as  opposed 
to  480  meter  spacing  in  figure  2. 10).  A  series  of  Stoneley  arrivals  appear  after  the  direct 
wave  and  have  a  phase  velocity  of  approximately  1.3  km/sec  (fig.  2.1  la).  Other  Stoneley 
waves  are  caused  by  scattering  from  the  P-diving  wave  (2.1  lb). 
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Figure  2.12.  Retrograde  elliptical  particle  motion  of  the  Stoneley  wave  shown  in  figure 
2. 1 1  b  at  7.5  km  range  for  a  seafloor  receiver. 
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Figure  2. 13.  Numerical  Schlieren  diagrams  for  the  other  models  with  Gaussian 
autocorrelation  functions  (see  figure  2.6  for  velocity  fields).  Only  the  4.0  second  snapshot 
is  shown  since  the  general  trends  and  principle  waves  are  similar  to  those  of  the  laterally 
homogeneous  case  (figure  2.7)  and  to  model  GAUSS50. 
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Figure  2.13.  Numerical  Schlieren  diagrams  for  the  other  models  with  Gaussian 
autocorrelation  functions  (see  figure  2.6  for  velocity  fields).  Only  the  4.0  second  snapshot 
is  shown  since  the  general  trends  and  principle  waves  are  similar  to  those  of  the  laterally 
homogeneous  case  (figure  2.7)  and  to  model  GAUSS50. 
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Figure  2.14.  Pressure  time  series  for  models  GAUSS  10  (2.14a),  GAUSS  100  (2.14b), 
and  GAUSS200  (2. 14c).  In  general,  the  amount  of  random  scatter,  or  'noise',  in  the 
seismograms  decreases  as  ka  increases.  Conversely,  deterministic  scattering  from  the 
random  heterogeneities  increases  as  ka  increases.  This  is  seen  by  an  increase  in  arrival  time 
and  amplitude  anomalies  as  well  as  the  appearance  of  unexpected  coherent  arrivals  (those 
not  seen  in  the  laterally  homogeneous  case)  in  the  seismograms. 
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Figure  2. 15.  a.  The  4.0  second  snapshot  for  model  SELFSIM200  with  a  self-similar 
velocity  autocorrelation  function.  This  model  contains  characteristics  of  the  Gaussian 
models  with  both  small  and  large  correlation  lengths.  The  small  scale  roughness  causes 
random  backscatter  and  the  generation  of  a  number  of  Stoneley  waves.  Larger  scale 
roughness  is  the  cause  of  coherent  secondary  body  wave  scattering  as  well  a  deterministic 
Stoneley  wave  scattering. 
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Figure  2. 15  b.  Pressure  time  series  for  model  SELFSIM200.  Again,  as  in  the  wavefront 
diagram,  this  model  contains  characteristics  of  scattering  from  both  large  and  small  scale 
heterogeneities. 
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DISCUSSION 


The  synthetic  seismograms  from  the  laterally  homogeneous  FLAT  model  (figure 
2.8)  contain  only  the  well  studied  deterministic  waves.  Scattering  from  lateral 
heterogeneities  causes  secondary  waves  to  complicate  the  seismic  traces  (see  figures  2.10, 
2.14,  and  2.15b).  For  ka  near  unity  these  secondary  arrivals  appear  as  random  incoherent 
'noise'  following  the  principal  arrivals.  Model  GAUSS50  (figure  2.10),  with  ka 
approximately  equal  to  one  at  the  seafloor,  contains  the  largest  amount  of  random  noise  for 
the  cases  presented  in  this  study.  It  has  been  shown  that  this  random  'noise'  is  caused  by 
the  scattering  of  body  waves  (both  P  and  S  types)  and  secondary  Stoneley  wave  generation 
from  interaction  of  the  principal  arrivals  with  the  heterogeneities.  An  inspection  of  the 
pressure  seismograms  in  figures  2.10  and  2.14  shows  that  in  general,  for  the  Gaussian 
models,  the  amount  of  random  'noise'  increases  as  ka  approaches  one. 

Deterministic  scattering  from  the  heterogeneities  begins  to  dominate  the  secondary 
features  of  the  seismograms  as  the  lateral  extent  of  the  heterogeneities  increases.  Model 
GAUSS  100  shows  much  less  of  the  random  'noise'  seen  in  GAUSS  10  and  GAUSS50. 
Instead,  more  coherent  secondary  arrivals  are  present  in  the  seismograms  of  GAUSS  100 
and  GAUSS200.  The  amount  of  incoherent  secondary  wave  energy  decreases  in  both 
GAUSS  100  and  GAUSS200  from  the  models  with  lower  ka  values.  In  fact,  almost  all  of 
the  secondary  arrivals  seen  in  GAUSS200  can  be  identified  as  coherent  secondary  Stoneley 
waves  or  secondary  scattered  body  waves  (mainly  secondary  S  waves  arriving  after  the 
principal  P-diving  wave  and  before  the  principal  S-diving  wave).  Larger  scale 
heterogeneities  also  have  a  greater  influence  on  the  amplitude  and  travel  time  anomalies  of 
the  principal  arrivals. 
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One  measure  of  the  influence  of  lateral  heterogeneities  on  the  coherence  of  wave 
propagation  is  the  cross  correlation  of  seismogram  arrivals  at  different  ranges.  Frankel  and 
Clayton  (1986)  carried  out  this  type  of  analysis  on  wave  propagation  through  a  halfspace 
with  random  heterogeneities.  Our  study  is  more  complicated  because  of  the  presence  of 
different  phases  arising  from  interaction  of  seismic  energy  with  the  liquid-solid  interface. 
Therefore,  before  we  can  compare  the  cross-correlation  patterns  of  different  models  with 
lateral  heterogeneity,  we  must  first  understand  the  behavior  of  the  cross-correlation  analysis 
within  the  laterally  homogeneous  FLAT  model  at  different  ranges  (and  thus  different  wave 
relationships). 

The  normalized  cross-correlation  coefficient  between  a  given  trace  and  any  other 
trace  is  equal  to  one  for  wave  propagation  through  a  homogeneous  box.  However, 
because  of  wave  interaction  with  the  liquid-solid  interface,  the  synthetic  seismograms  in 
this  study  contain  many  different  wave  types  on  a  given  seismic  trace.  Because  the 
different  waves  travel  with  different  speeds  and  directions  and  the  source  is  not  equidistant 
from  all  of  the  receivers,  a  cross-correlation  between  any  two  entire  traces  is  going  to 
contain  these  deterministic  effects.  As  the  distance  between  traces  increases,  the  cross- 
correlation  coefficient  should  roughly  decrease  because  of  unique  characteristics  of  the 
different  waves  present.  We  have  chosen  to  look  at  only  the  first  0.5  seconds  of  the  initial 
P-wave  arrival  (the  direct  wave  at  small  ranges,  the  P-diving  wave  at  large  ranges)  for  the 
cross-correlation  analysis  in  order  to  minimize  the  artifacts  of  differences  in  the  wave  types. 

Figure  2.16  is  a  plot  of  the  cross-correlation  coefficients  between  a  number  of 
traces  (0.5  seconds  after  the  P-diving  wave  first  zero  crossing)  which  are  the  seafloor 
pressure  time  series  of  model  FLAT  from  6.0  and  9.0  km.  (see  labels  on  curves  in  figure 
2.16)  and  all  of  the  other  traces  within  this  synthetic  seismogram  between  6.0  and  9.0  km. 
This  figure  demonstrates  that  the  correlation  between  P-diving  waveforms  is  range 
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Figure  2.16.  Cross  correlation  coeficients  for  the  synthetic  seismogram  traces  between  6.0 
and  9.0  km.  of  model  FLAT.  The  curves  shown  are  for  a  number  of  reference  traces 
(labelled)  correlated  with  all  other  traces  between  6.0  and  9.0  km.  Spatial  lag  is  the 
difference  between  the  range  of  the  reference  trace  and  that  of  the  cross-correlation  trace. 
Correlation  analysis  was  done  for  the  first  0.5  seconds  after  the  P-diving  wave  first  zero 
crossing.  Note  that  the  correlation  coefficient  falls  off  more  rapidly  for  reference  traces  at 
relatively  greater  ranges. 
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Although  significant  'coda'  exists  for  model  GAUSS  10,  suggesting  that  random 
scatter  is  important,  there  is  apparendy  little  change  in  the  waveform  caused  by  the  small 
scatterers.  Since  ka  is  small  for  GAUSS  10,  the  wave  field  is  starting  to  see  the  velocity 
perturbations  as  an  equivalent  body  with  an  average  velocity  field  approaching  that  of  the 
FLAT  model.  The  random  noise  is  slightly  higher  in  frequencey  and  arises  from  scattering 
of  the  higher  frequency  components  of  the  source.  The  predominantly  10  Hz  P-diving 
wave  is  not  greatly  affected. 

The  distinction  between  the  cross  correlation  characteristics  of  the  different  models 
decreases  as  the  range  of  the  reference  trace  increases.  A  clear  separation  exists  between 
the  cross-correlation  patterns  of  the  different  models  if  the  reference  trace  is  at  a  relatively 
small  range  (such  as  6.0  km.  in  figure  2.17a.).  As  the  range  of  the  reference  trace 
increases  there  is  an  overall  decrease  in  most  of  the  absolute  values  of  the  correlation 
coeficients  as  well  as  a  decrease  in  the  distinction  of  the  patterns  between  models.  This  is 
probably  due  to  the  influence  of  the  interference  head  wave  at  greater  ranges.  That  is,  low 
correlation  between  traces  is  influenced  more  by  the  interference  head  wave  than  by  the 
heterogeneity  at  large  ranges.  Model  GAUSS50  which  has  a  lower  correlation  curve  than 
the  others  has  about  the  same  correlation  values  with  different  reference  ranges  because  the 
large  amount  of  scattering  causes  the  interference  head  wave  to  be  at  a  low  amplitude 
already  (see  figures  2. 10  and  2. 14  for  a  comparison  of  interference  head  wave  amplitudes). 
Models  with  less  random  scatter  have  stronger  interference  head  waves  at  large  ranges  so 
that  the  decrease  in  correlation  coefficients  with  range  is  larger  for  these  models. 


163 


Figure  2. 17.  Comparison  of  the  cross-correlation  curves  for  all  Gaussian  models  and  the 
laterally  homogeneous  model  at  4  different  reference  trace  ranges,  a.  Reference  trace  at 
6.0  km.  There  is  a  clear  distinction  between  the  curves  for  the  different  models  at  this 
range,  b.  7.0  km.  reference  trace,  c.  8.0  km.  reference  trace,  d.  9.0  km.  reference 
treace.  The  absolute  values  for  the  GAUSS50  curve  have  remained  about  the  same  for  the 
4  different  reference  ranges  presented.  However,  the  overall  values  for  the  other  curves 
have  decreased  with  increasing  range  so  that  by  9.0  km.  it  is  difficult  to  notice  any  major 
distinctions  between  the  5  curves. 
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CONCLUSIONS 


We  have  demonstrated  the  applicability  of  the  finite  difference  method  to  forward 
modeling  of  elastic  wave  scattering  from  random  lateral  heterogeneities  in  the  upper  oceanic 
crust.  For  random  velocity  perturbations  with  Gaussian  autocorrelation  functions  there  is  a 
general  trend  of  increasing  random  scatter  in  the  seismograms  as  ka  approaches  one. 
Secondary  deterministic  scattering  from  the  heterogeneities  increases  as  ka  increases  from 
one.  Velocity  perturbations  with  a  self-similar  autocorrelation  function  cause  scattering 
with  characteristics  of  both  the  large  and  small  correlation  length  Gaussian  models.  That 
is,  the  self-similar  model  contains  the  random  backscatter  seen  in  the  Gaussian  models  with 
ka  less  than  one  as  well  as  the  effects  of  travel  time  and  amplitude  anomalies  seen  in  the 
Gaussian  models  with  ka  larger  than  one. 

Heterogeneities  in  the  upper  crust  near  the  water-solid  interface  act  as  sources  for 
secondary  Stoneley  waves  when  illuminated  by  primary  seismic  energy.  The  presence  of 
the  water-solid  interface  in  the  models  is  important  because  the  majority  of  the  seafloor 
'noise'  seen  in  the  synthetic  seismograms  seems  to  travel  as  secondary  Stoneley  waves. 
Simple  models  of  propagation  through  a  half-space  (e.g.  Frankel  and  Clayton,  1986, 
McLaughlin  et.al.,  1985)  do  not  include  interface  phenomena  which  are  an  integral  part  of 
the  marine  problem.  Also,  inclusion  of  the  water-solid  interface  allows  the  examination  of 
the  influence  of  scattering  on  the  many  different  wave  types  which  arise  from  interaction 
with  the  interface.  The  interference  P-head  wave  has  a  significant  influence  on  the  cross¬ 
correlation  between  traces  in  the  scattering  models. 

Cross-correlation  analysis  of  the  P-diving  wave  in  the  scattering  models  reveals  two 
important  characteristics.  First,  as  ka  approaches  one  from  below,  the  correlation  between 


166 


seismic  traces  decreases.  Deterministic  scattering  from  spatially  larger  heterogeneities 
ikii>  \ )  causes  much  less  decorrelation  of  the  traces  than  does  'random'  scattering  from 
heterogeneities  with  ka  near  one.  Second,  the  cross-correlation  values  are  dependent  upon 
the  range  of  the  reference  trace  of  the  analysis.  This  is  due  to  the  influence  of  the 
interference  head  wave  at  large  ranges. 
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Seismo/acoustic  propagation  through  rough  seafloors 


INTRODUCTION 


The  problem  of  seismic  energy  scattering  from  heterogeneities  within  the  cmst  is 
complicated  by  the  presence  of  topography  on  the  seafloor  or  sediment/basement  interface 
While  heterogeneities  within  the  crust  can  contain  velocity  contrasts  on  the  order  of  1-20 
per  cent,  the  water-solid  interface  can  represent  a  velocity  contrast  of  over  100  per  cent. 
Because  of  this  sharp  impedance  contrast,  non-planar  seafloors  can  cause  large  scattering 
effects  when  energy  propagates  into  the  oceanic  crust.  With  increased  interest  in  the 
seismic  study  of  mid-ocean  ridges  where  the  impedance  contrast  is  particularly  high,  it  is 
important  to  understand  the  relationship  between  topography  and  seismo/acoustic 
propagation  through  the  crust. 

Topography  of  the  seafloor  covers  a  wide  spectrum  of  wavenumbers.  While  the 
general  impression  of  the  seafloor  is  one  of  a  deep,  flat,  sediment  covered  abyssal  plane, 
Menard  (1964)  points  out  that  most  of  the  world  ocean  floor  is  indeed  covered  by  rough, 
hilly  terrain  (as  much  as  85%  of  the  Pacific  ocean  floor  may  have  this  type  of  terrain). 

Early  surveys  of  compilations  of  the  statistics  of  ocean  floor  topography  showed  a  very  red 
spectrum  (Bell,  1975;  Bell,  1979),  that  is,  one  which  is  weighted  heavily  in  the  longer 
wavelengths.  However,  inclusion  of  more  recent  deep-tow  and  mid-ocean  ridge  data  will 
add  at  least  some  power  to  the  spectrum  at  higher  wavenumbers.  Regardless,  there  is 
enough  topography  with  wavenumbers  similar  to  those  of  seismic  energy  to  merit  study  of 
the  propagation  effects  of  topography.  No  section  of  the  ocean  floor  is  perfectly  smooth, 
but  rather,  each  contains  a  range  of  topographic  length  scales.  The  section  of  the  range  of 
length  scales  seen  will  depend  on  the  resolution  of  the  detection  method  used  (Akal,  1978; 
Ogilvy,  1987).  Seismic  energy  will  be  most  effective  in  detecting  topography  which  has 
length  scales  similar  to  and  larger  than  that  of  the  seismic  wavelength. 
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Acoustic  scattering  from  rough  surfaces  is,  in  general,  a  very  well  studied  problem. 
Many  approximations  exist  to  calculate  the  field  scattered  from  a  rough  bottom  (Ogilvy, 
1987).  Well  known  solutions  which  contain  restrictive  approximations  include  Rayleigh’s 
method  of  plane  wave  summation  (Rayleigh,  1878),  the  method  of  small  perturbations 
(very  small  topography  with  respect  to  acoustic  wavelength),  the  Kirchhoff  method 
(assumes  shallow  slopes),  and  methods  which  sum  contributions  from  distributions  of 
points  sources,  facets,  or  protuberances  along  the  bottom  (see  Ogilvy,  1987). 

Recently,  methods  have  been  developed  which  are  more  flexible  and  require  fewer 
limiting  assumptions  and  give  more  realistic  results.  Bouchon  (1985)  used  a  boundary 
integral  formulation  coupled  with  the  discrete  wavenumber  method  to  calculate  the  scattered 
field  emanating  from  a  rough  surface  of  arbitrary  wavenumber.  This  method  has  been 
successfully  used  to  model  diffraction  of  energy  from  a  corrugated  boundary  (Campillo  and 
Bouchon,  1985;  Paul  and  Campillo,  1988).  The  emergence  of  more  powerful  computers 
has  made  the  use  of  numerical  methods  such  as  the  finite  difference  and  finite  element 
methods  much  more  attractive  because  of  the  complete  solution  to  the  scattering  problem 
with  very  few  limiting  assumptions.  The  major  limitation  of  these  methods  is  that  they  are 
computationally  intensive  and  unpractical  for  very  long  range  models.  However,  these 
methods  are  well  suited  for  investigating  the  local  effects  of  scattering  at  the  sea  floor. 
Stephen  (1984)  used  the  finite  difference  method  to  explain  double  head  waves  caused  by 
seafloor  topography  and  Dougherty  and  Stephen  (1987,  chapter  1)  used  this  method  to 
study  'refraction  branch  diffractions'  from  seafloor  features  beyond  the  critical  range.  (Hill 
and  Levander,  1984;  Levander  and  Hill,  1985)  also  used  the  finite  difference  method  and 
investigated  scattering  from  rough  low  velocity  layers  within  the  crust  and  rough  elastic 
surface  layers. 
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The  work  for  this  chapter  involves  the  use  of  a  velocity-stress  formulation  of  the 
finite  difference  method  for  elastic  wave  propagation  (Dougherty  and  Stephen,  1988 
(chapter  2);  Stephen,  1988;  Virieux,  1986)  to  investigate  scattering  from  rough  seafloor 
topography.  While  the  long  term  goal  of  this  work  is  to  be  able  to  distinguish  between 
scattering  from  volume  and  surface  heterogeneities,  this  chapter  deals  mainly  with 
procedural  problems  with  the  method.  Briefly,  most  of  the  difficulty  in  using  this  method 
for  modelling  smoothly  varying  topography  lies  in  the  definition  of  a  non-planar, 
curvilinear  seafloor  on  a  rectangular  grid.  For  rectangular  grids  which  are  numerically 
stable,  the  discretization  of  the  seafloor  profile  imposes  a  stair-stepping  microroughness 
onto  the  larger  scale  sinusoidal  topography.  Scattering  phenomena  due  to  this  problem  will 
be  discussed  as  well  as  'real'  effects  which  can  be  expected  from  sinusoidal  topography 
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INITIAL  MODELS 


The  goal  of  this  and  future  work  on  scattering  from  rough  surfaces  is  to  provide 
insight  into  the  distinction  between  topographic  and  volumetric  scattering  mechanisms  and 
effects.  There  are  obviously  an  infinite  number  of  topographic  profiles  present  in  the 
oceans  so  we  have  no  hope  of  modelling  every  specific  case.  The  best  that  can  be  done  is 
to  quantify  characteristics  of  scattering  from  generalized  seafloor  structures.  In  an  effort  to 
start  this  process,  this  chapter  deals  with  a  number  of  rough  seafloor  models  which  contain 
sinusoidal  topography  with  height  on  the  order  of  the  seismic  wavelength.  While  the 
existence  of  exactly  this  type  of  topography  is  highly  unlikely,  these  models  do  provide  a 
good  starting  point  for  the  problem  and  are  useful  in  that  they  point  out  both  some 
important  scattering  mechanisms  as  well  as  operational  constraints  on  the  method  used. 

Sinusoidal  seafloor  profiles  were  chosen  for  this  first  set  of  models  for  a  number  of 
reasons.  The  regular  nature  of  the  sine  wave  makes  the  topography  and  underlying 
velocity  structure  easy  to  define  on  the  finite  difference  grid.  Also,  the  regularity  of  the 
sine  functions  allows  for  easy  modification  of  height  and  wavelength  parameters.  The 
disadvantage  of  using  a  regular  structure  is  that  Bragg  interference  effects  may  become  an 
issue.  Another  reason  why  the  sine  wave  was  chosen  was  because  it  provides  for  an  easily 
defined,  smoothly  varying  bottom.  The  effects  of  sharp  breaks  (such  as  the  quarter  space 
problem)  have  been  investigated  by  a  number  of  different  authors  (Achenbach  et  al„  1982; 
Kelly  et  al.,  1976;  Virieux,  1984;  Virieux,  1986).  Finally,  in  a  very  rough  sense, 
sinusoidal  topography  does  approximate  linear,  ridge-parallel,  seafloor  structures  near  mid¬ 
ocean  ridges. 
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Elastic  wave  propagation  through  a  number  of  rough  seafloor  models  was 
calculated  by  using  a  2-D  elastic  finite  difference  formulation  of  the  elastodynamic  wave 
equations  for  heterogeneous  media.  This  formulation  is  well  suited  to  handle  laterally 
heterogeneous  media  and  the  interface  between  a  liquid  and  solid  present  in  marine  models 
(Dougherty  and  Stephen,  1988  (chapter  2);  Stephen,  1988b;  Virieux,  1986).  More  detail 
on  the  specific  formulation  used  can  be  found  in  Dougherty  and  Stephen  (1988,  chapter  2) 
and  Stephen  (1988a).  A  diagram  of  ihe  grid  layout  is  presented  in  figure  3.1.  A 
rectangular  grid  with  vertical  and  norizontal  spacing  of  10  meters  was  used  for  most  of  the 
models  presented  in  this  chapter.  Outgoing  energy  is  absorbed  on  the  top,  bottom  and  right 
hand  sides  of  the  grid  (Dougherty  and  Stephen,  1988  (chapter  2);  Stephen,  1988c,  Clayton 
and  Engquist,  1977),  while  the  left  hand  side  acts  as  an  axis  of  symmetry.  The  sinusoidal 
profiles  have  a  phase  shift  of  -rt/2  to  avoid  symmetry  problems  at  zero  range.  Two  rows  of 
receivers  (hydrophones)  are  located  20  meters  above  the  seafloor  and  500  meters  above  the 
topography  baseline  (see  figure  3.1).  The  row  of  receivers  located  just  above  the  seafloor 
is  not  horizontal,  but  rather,  it  follows  the  topography  of  the  seafloor. 

Initially,  a  set  of  sinusoidal  models  with  10  meter  grid  spacing  and  a  variety  of 
amplitude  and  wavelength  values  were  run.  These  seafloor  profiles  are  shown  in  figure 
3.2  and  the  velocity-depth  function  of  model  SIN360  for  average,  hill,  and  valley  sections 
are  shown  in  figure  3.3.  The  'average'  velocity-depth  function  used  is  the  same  as  for  the 
Dougherty  and  Stephen  (1987,  1988,  chapters  1  and  2)  FLAT  model  and  is  typical  of 
young  ocean  crust  (Ewing  and  Purdy,  1982;  Purdy,  1982a).  The  effect  of  hills  and  valleys 
on  the  velocity  profiles  is  to  decrease  the  second  gradient  under  hills  and  increase  the 
second  gradient  under  valleys. 

The  two-dimensional  P-wave  velocity  profile  for  model  SIN360  is  shown  in  figure 
3.4.  This  model  has  sinusoidal  seafloor  topography  with  a  wavelength  of  3.6  kilometers 
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and  amplitude  of  150  meters  (300  meters  from  peak  to  valley).  Features  of  the  results  of 
this  model  demonstrate  both  the  characteristics  of  scattering  from  large  deterministic 
seafloor  structures  as  well  as  the  limitations  of  the  finite  difference  grid  (with  coarse  grid 
spacing)  in  handling  the  non-planar  water-solid  interface. 

As  seen  in  figures  3.3  and  3.4,  the  gradient  in  P-wave  velocity  is  different  between 
hills  and  valleys.  In  general,  the  first  segment  of  the  'average'  velocity-depth  function 
follows  the  seafloor,  that  is,  the  vertical  gradient  (not  the  normal  gradient)  is  the  same  at  all 
points  in  range.  The  second  gradient  segment  is  deformed  under  the  topography  to  be 
larger  under  valleys  and  smaller  under  hills.  Also,  the  effective  gradient  for  waves 
normally  incident  on  the  topography  will  be  greater  than  the  actual  vertical  gradient.  Both 
of  these  factors  will  add  to  the  travel  time  anomalies  present  in  the  time  series. 

Angle  of  incidence  for  a  P-wave  source  at  zero  range  is  plotted  against  range  at  the 
seafloor  for  model  SIN360  in  figure  3.5.  Also  shown  on  the  figure  are  the  critical  angles 
and  ranges  for  P  and  S-diving  waves  from  the  surface  source.  Critical  ranges  for  P  and  S- 
diving  waves  in  figure  3.5  occurs  where  the  angle  of  incidence  curve  crosses  the  P  or  S- 
diving  wave  critical  angle  line.  Backscattering  will  occur  when  angle  of  incidence  is  zero 
and  ray  theory  'shadow  zones'  would  occur  beyond  ranges  where  the  source  path  is 
tangent  to  the  seafloor  (angle  i  =  tt/2).  Of  course,  no  actual  'shadow  zones'  are  rt  -esent 
because  of  'Franz-type'  waves  which  are  diffracted  from  the  curved  surface  at  grazing 
incidence.  These  diffracted  waves  are  analogous  to  the  circumferential  waves  seen  in 
studies  of  acoustic  plane  wave  scattering  from  large  elastic  cylinders  (Neubauer,  1986). 
Multiple  head  waves  for  both  compressional  and  shear  waves  will  occur  whenever  the 
angle  of  incidence  curve  crosses  the  critical  angle  curve  for  the  respective  wave  types  then 
falls  back  below  it  (Stephen,  1984). 
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Figure  3.1.  Finite  difference  grid  parameters  for  the  initial  rough  seafloor  models.  The 
seafloor  is  located  at  a  mean  depth  of  3.0  km.  (topography  baseline).  Hydrophones  are 
located  every  40  meters  in  range  along  one  horizontal  row  in  the  water  column  and  one  row 
which  follows  the  seafloor  20  meters  above  the  topography. 
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Figure  3.2.  Sinusoidal  seafloor  topography  profiles  for  the  four  initial  models. 
Topography  wavelengths  for  the  models  are  900  meters  (a.,  b.),  1.8  km.  (c.),  and  3.6  km. 
(d.)  with  amplitudes  of  80  meters  (a.)  and  150  meters  (b.,  c.,  d.). 
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Figure  3.3.  Velocity-depth  functions  for  model  SIN360  beneath  hills  (short  dashes), 
valleys  (long  dashes)  and  average  structure  (solid  line).  The  first  gradient  segment  is 
constant  regardless  of  structure  and  the  second  gradient  segment  is  either  lengthened  or 
compressed  to  meet  a  constant  velocity  of  6.0  km/sec  at  4.3  km.  depth.  Shear  wave 
velocity-depth  profiles  are  similar  and  assume  a  Poisson's  ratio  of  0.25. 
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Figure  3.4.  Compression^  wave  velocity  profile  for  model  SIN360.  Horizontal  gradients 
are  imposed  by  the  topography,  as  well  as  vertical  gradients  which  vary  with  range.  The 
shear  wave  velocity  profile  is  similar  and  assumes  a  Poisson's  ratio  of  0.25. 
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Figure  3.5.  Angle  of  incidence  vs.  range  for  a  pulse  source  at  0.0  km.  range  and  3.0  km. 
above  the  seafloor  for  model  SIN360.  Topography  allows  for  multiple  S-wave  critical 
ranges.  Theoretically,  every  range  where  the  critical  angle  falls  below  the  critical  angle  for 
P  or  S  energy  (given  by  the  dashed  lines),  P  or  S  energy  can  be  directly  transmitted  into  the 
bottom. 
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Model  SIN360,  as  well  as  models  SIN90A,  SIN90B,  and  SIN  100  presented 
below,  were  run  for  5000  time  steps  at  one  millisecond  per  time  step  for  a  total  real  time  of 
5  seconds  after  the  surface  shot  (and  a  10  Hz  source).  Numerical  Schlieren  diagrams 
(figure  3.6)  were  output  every  400  time  steps  (0.4  seconds)  starting  at  2.0  seconds.  These 
diagrams  represent  normalized  compressional  (top)  and  shear  (bottom)  transient  energy 
within  the  model  environment  (chapter  2  for  more  on  this  format).  Dougherty  and  Stephen 
(1988,  chapter  two)  present  both  the  snapshots  and  time  series  in  the  same  format  as  this 
study  for  the  laterally  homogeneous  model  (FLAT)  with  a  horizontal,  planar  seafloor.  This 
laterally  homogeneous  FLAT  model  can  be  considered  as  a  'control'  model  for  the  work 
shown  here. 

Figure  3.6  is  a  series  of  snapshots  which  represents  the  time  evolution  of 
wavefronts  (P-waves  top  and  shear  waves  bottom)  within  model  SIN360.  The  initial  P- 
wave  pulse  source  is  just  above  the  seafloor  after  2.0  seconds  (figure  3.6a).  At  the  2.4 
second  snapshot  (figure  3.6b)  the  wavefront  has  partitioned  into  transmitted  P  and  S, 
reflected  P,  and  direct  P  waves.  Note  particularly  the  curved  P-wave  front  which  has  been 
reflected  from  the  left  side  of  the  first  hill.  This  will  also  show  up  as  additional  curvature  in 
the  P-diving  wave  arrival  in  the  t-x  space  of  the  time  series.  By  2.8  seconds  after  the  shot, 
the  wave  interactions  have  gotten  much  more  complex  (figure  3.6c).  Multiple  shear  and 
compressional  transmitted  waves  at  very  close  range  are  probably  caused  by  constructive 
interference  between  wavefronts  scattered  from  the  stepwise  curved  seafloor  (the  curved 
interface  is  made  up  of  a  number  of  10  meter  rectangular  steps).  It  will  be  shown  later  that 
when  the  grid  increment  is  made  smaller,  much  of  this  scattering  will  disappear. 

There  are  a  number  of  notable  scattering  features  which  appear  in  the  snapshots 
from  3.2  to  4.8  seconds  (Figures  3.6d-h).  Most  prominent  among  these  is  the 
preponderance  of  scattered  energy  which  is  found  between  the  primary  arrivals  (those  seen 
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Figure  3.6  a,b.  Snapshots  of  normalized  compressional  (top)  and  shear  (bottom)  energy 
for  model  SEN360.  The  sinusoidal  seafloor  is  located  by  the  solid  line  at  an  average  depth 
of  3.0  km.  3.6a.  2.0  seconds  after  the  shot.  The  cylindrical  wave  from  the  pulse  line 
source  at  the  sea  surface  is  just  reaching  the  seafloor  at  zero  range.  No  shear  energy  is 
present  yet  in  the  system.  3.6b.  2.4  seconds  after  the  shot.  Seismic  energy  has  interacted 
with  the  seafloor  and  partitioned  into  direct,  seafloor  reflected,  and  transmitted  waves  of 
both  compressional  (top  frame)  and  shear  (bottom  frame)  types. 
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Figure  3.6  c,d.  Snapshots  of  normalized  compressional  (top)  and  shear  (bottom)  energy 
for  model  SIN360.  The  sinusoidal  seafloor  is  located  by  the  solid  line  at  an  average  depth 
of  3.0  km.  Times  of  2.8  seconds  (3.6c.)  and  3.2  seconds  (3.6d.)  after  the  surface  shot. 
Shear  energy  will  continue  to  be  transmitted  into  the  crust  out  to  a  range  of  approximately 
5.5  kilometers  (see  Figure  3.5).  The  energy  which  appears  random  and  incoherent  has 
been  scattered  from  the  small  grid-induced  steps  of  the  sloping  topography. 
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Figure  3.6  e,f.  Snapshots  of  normalized  compressional  (top)  and  shear  (bottom)  energy 
for  model  SIN360.  The  sinusoidal  seafloor  is  located  by  the  solid  line  at  an  average  depth 
of  3.0  km.  Times  of  3.6  seconds  (3.6e.)  and  4.0  seconds  (3.6f.)  after  the  surface  shot. 
The  curved  arrivals  in  the  shear  plot  of  4.0  second  (3.6f.)  around  range  2.0  km.  are 
travelling  to  the  left  and  are  shear  waves  which  have  been  backscattered  from  the  rough 
seafloor.  These  show  up  as  curved  arrivals  on  VSP  receivers. 
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Figure  3.6  g,h.  Snapshots  of  normalized  compressional  (top)  and  shear  (bottom)  energy 
for  model  SIN360.  The  sinusoidal  seafloor  is  located  by  the  solid  line  at  an  average  depth 
of  3.0  km.  Times  of  4.4  seconds  (3.6g.)  and  4.8  seconds  (3.6h.)  after  the  surface  shot. 
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There  are  a  number  of  notable  scattering  features  which  appear  in  the  snapshots 
from  3.2  to  4.8  seconds  (figures  3.6d-h).  Most  prominent  among  these  is  the 
preponderance  of  scattered  energy  which  is  found  between  the  primary  arrivals  (those  seen 
in  the  FLAT  model  of  chapters  one  and  two  plus  additional  diving  waves  associated  with 
the  sloping  sea  floor).  Since  the  slopes  in  model  SIN360  are  relatively  small,  small  scale 
diffraction  type  scattering  is  not  expected.  Therefore,  the  incoherently  scattered  energy 
must  be  due  to  the  steps  imposed  by  coarse  sampling  of  the  topography.  This  energy  is 
scattered  in  all  directions  from  the  seafloor.  Also,  a  significant  amount  of  energy  is 
backscattered  or  reverberated  into  the  water  column.  The  accuracy  of  this  scattering  will  be 
discussed  in  more  detail  below,  but  .'or  these  preliminary  models,  it  appears  that  small 
(relative  to  the  seismic  wavelength)  sharp  topographic  features  can  have  a  significant  effect 
on  the  scattering  of  seismo/acoustic  energy.  One  consequence  of  this  scattering  at  the 
stepwise  interface  is  that  interface  waves  can  be  generated  at  the  boundary.  These  waves 
have  the  characteristics  of  Stoneley  waves  and  appear  as  small  packets  of  energy  confined 
to  the  immediate  area  of  the  seafloor  (amplitude  decays  exponentially  away  from  both  sides 
of  the  seafloor)  and  they  travel  at  a  velocity  just  less  than  the  compressional  wave  velocity 
in  the  water  (Luppe  and  Doucet,  1988;  Scholte,  1947;  Stoneley,  1924).  These  interface 
waves  are  similar  to  those  generated  by  scattering  from  volume  heterogeneities  with  a  Hat 
seafloor. 

Compressional  and  shear  diving  waves  generated  from  an  acoustic  pulse  source  in  a 
planar  seafloor  environment  are  singular  and  continuous.  However,  when  the  topography 
on  the  seafloor  is  steep  enough  to  cause  the  angle  of  incidence  ot  the  direct  wave  to 
repeatedly  fall  below  the  critical  angle  for  diving  energy,  additional  distinct  P  and  S-diving 
waves  can  occur.  These  additional  diving  waves  must  not  be  confused  with 
compressional  and  shear  diving  wave  multiples  which  are  caused  by  primary  diving  wave 
bounces  off  the  seafloor  from  below.  This  critical  angle  induced  process  may  be  an 


194 


of  incidence  to  fall  below  the  critical  angle  for  P-waves  after  the  primary  P-wave  (see  figure 
3.5).  There  is.  however,  a  relatively  strong  P-wave  multiple  that  can  be  seen  in  the  3.2 
second  snapshot  (figure  3.6d)  which,  along  with  the  interference  P-headwave  (made  up  of 
a  superposition  of  all  P-wave  multiples)  rapidly  decreases  in  amplitude  in  succeeding 
snapshots.  The  shear  wave  critical  angle  is  much  higher  than  that  for  compressional  waves 
and  the  topography  of  SIN360  is  steep  enough  to  allow  for  an  additional  S-diving  wave. 
This  diving  wave  falls  directly  behind  the  primary  S-diving  wave  in  the  4.0  second 
snapshot. 

A  more  traditional  means  of  observing  the  synthetic  seismic  data  is  by  the  use  of 
time  series  of  hydrophone  response.  Hydrophone  receivers  were  placed  in  2  rows,  one 
located  20  meters  above  the  seafloor  and  one  horizontal  array  located  500  meters  above  the 
topography  baseline  (the  point  in  the  sine  curve  where  amplitude  is  zero)  as  seen  in  figure 
3.1.  Hydrophones  were  spaced  every  8  grid  points  or  80  meters  apart  in  both  receiver 
rows.  The  time  series  from  both  rows  of  receivers  of  model  SEN360  are  shown  in  figure 
3.7.  Two  amplitude  scales  are  presented  for  each  seismogram.  The  smaller  amplitude 
scale  was  used  to  display  an  unclipped  direct/reflected  wave  and  a  scale  factor  five  times 
larger  was  used  to  display  a  visible  and  unclipped  P-diving  wave.  None  of  the  seismogram 
arrivals  have  been  corrected  for  topography  (as  in  Purdy,  1982b)  so  that  in  general  the 
arrivals  (except  for  the  direct  arrival)  for  the  horizontal  row  of  receivers  will  appear  more 
curved  than  those  for  the  seafloor  hydrophones. 

Figures  3.7c  and  3.7d  are  the  pressure  time  series  for  the  row  of  hydrophones 
along  the  seafloor  plotted  at  the  two  different  amplitude  scales.  The  P-diving  wave  is  only 
slightly  curved  due  to  travel  time  anomalies  induced  by  the  different  velocity-depth  profiles 
below  hills  and  valleys.  Amplitude  anomalies  due  to  focussing/defocussing  around  the 
topography  are  prominent  along  the  entire  range  of  the  arrival  when  compared  to  the  FLAT 
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Figure  3.7  a,b.  Waterborne  hydrophone  receivers  for  model  SIN360.  3.7a.  Arrivals  have 
been  scaled  to  show  an  unclipped  direct  wave.  Discontinuity  of  the  reflected  arrival  is 
caused  by  multiple  shear  wave  critical  ranges  caused  by  the  sinusoidal  topography  (see 
figure  5).  3.7b.  Amplitudes  are  increased  by  a  scale  factor  of  five  over  figure  3.7a.  to 
show  the  unclipped  P-diving  wave.  Scattered  arrivals  after  the  reflected  and  P-diving 
waves  are  quite  significant  with  respect  to  the  P  and  S  primary  transmitted  waves. 
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Figure  3.7  c,d-  Seafloor  hydrophone  receivers  for  model  SIN360.  3.7c.  Amplitudes 
relative  to  unclipped  direct  wave.  Interface  waves  occur  following  the  direct  wave.  This 
'noise'  does  not  appear  on  the  waterborne  receivers.  3.7d.  Amplitudes  relative  to  the  P- 
diving  wave.  The  signal  following  the  direct  wave  is  swamped  by  scattered  interface  and 
body  waves.  Energy  scattered  out  of  the  P-diving  wave  by  the  topography  arrives  between 
the  P  and  S  primary  waves. 
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model  of  Dougherty  and  Stephen  (1987,  1988,  chapters  1  and  2).  As  discussed  earlier, 
there  are  no  additional  P-diving  waves  and  the  interference  P-wave  is  also  not  present.  The 
lack  of  occurrence  of  the  interference  P-wave  is  due  to  destructive  interference  caused  by 
the  rough  topography  (Dougherty  and  Stephen,  1987  (chaptet  1)).  Additional  shear 
waves,  on  the  other  hand,  do  appear  in  figure  3.7d.  In  fact,  the  second  shear-diving  wave 
has  a  stronger  amplitude  than  the  first  arrival  at  some  ranges.  Scattered  shear  waves  are 
present  both  after  the  direct  wave  (although  these  show  up  better  in  the  waterborne 
receivers)  and  between  the  P-diving  wave  and  the  S-diving  waves  (primary  and  additional 
arrivals). 

The  most  obvious  feature  of  the  seafloor  receivers  is  the  amount  of  strong  energy 
arriving  after  the  direct  wave  for  ranges  between  0  and  6  kilometers.  While  some  of  this 
energy  is  made  up  of  scattered  body  waves,  the  majority  of  the  energy  is  travelling  in  the 
form  of  an  interface  wave.  The  Stoneley  waves  are  travelling  in  both  the  forward  and 
reverse  directions  but  are  difficult  to  distinguish  from  the  direct  wave  in  the  forward 
direction.  Some  Stoneley  waves  are  also  generated  by  scattering  of  the  P-diving  wave  and 
show  up  just  later  in  time  than  the  P-diving  arrival  (i.e.  at  approximately  7.5  km.  in  figure 
3.7d).  As  is  expected,  the  Stoneley  waves  have  a  phase  velocity  less  than  the  acoustic 
velocity  in  the  water  and  in  general  have  very  large  amplitude,  comparable  in  some  places 
to  the  amplitude  of  the  direct  compressional  wave. 

Most  of  the  features  present  in  the  seafloor  seismograms  (figures  3.7c  and  3.7d)  are 
also  present  on  the  waterborne  receivers  seen  in  figures  3.7a  and  3.7b.  Again,  these 
receivers  are  in  a  horizontal  row  located  500  meters  above  the  topography  baseline.  The 
most  obvious  difference  between  seismograms  of  the  two  rows  of  receivers  is  the  relative 
lack  of  'noise'  or  energy  arriving  after  the  seafloor  reflection  of  the  direct  wave.  This 
energy  in  the  seafloor  time  series  is  the  Stoneley  wave  energy  which  decreases 
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exponentially  in  amplitude  away  from  the  seafloor.  The  row  of  receivers  placed  500  meters 
above  the  seafloor  does  not  pick  up  this  energy  confined  to  the  seafloor.  As  a  result,  the 
backscattered  body  waves,  mainly  of  S-type,  show  up  as  distinct  arrivals  after  the  reflected 
wave  in  the  upper  row  of  receivers.  In  general,  however,  the  amount  of  backscattered 
energy  is  much  less  in  these  receivers  than  for  those  along  the  seafloor.  The  greater 
curvature  of  the  P-diving  and  reflected  waves  is  due  to  variable  water  path  lengths  caused 
by  the  topography.  If  travel  times  were  explicitly  being  investigated,  this  could  be 
corrected  in  part  by  applying  a  topography  correction  to  the  data. 

Another,  way  to  observe  body  waves  within  the  crust  is  the  vertical  seismic  profile 
(VSP).  A  VSP  is  carried  out  by  recording  normally  incident  shots  at  a  number  of  locations 
in  a  borehole  or  along  a  vertical  array.  When  the  source  and  receiver  array  are  offset  from 
normal  incidence,  this  type  of  experiment  is  referred  to  as  an  oblique  seismic  experiment 
(OSE,  Stephen,  1983;  Stephen  et  al.,  1979;  Stephen.,  1980;  Stephen,  1984).  By  using 
the  finite  difference  method,  we  have  the  ability  to  place  receivers  anywhere  on  the  grid. 
SIN360  was  rerun  with  a  vertical  arrays  of  geophones  starting  at  zero  range  and  spaced 
every  kilometer  out  to  the  boundary  of  the  model.  Scattered  body  waves  will  have  a  similar 
appearance  in  both  the  VSP  and  OSE  geomcudes.  The  purpose  here  for  looking  at  a 
vertical  array  is  not  to  obtain  velocity  information,  but  rather,  to  observe  sideways 
travelling  waves  (specifically  backscattered  body  waves). 

The  vertical  array  of  geophones  can  be  seen  in  figure  3.1  and  a  cartoon  of 
'expected'  P-wave  arrivals  for  a  laterally  homogeneous  environment  is  shown  in  figure 
3.8.  A  reflected  arrival  from  a  layer  within  the  depth  of  the  vertical  array  would  intersect 
one  of  the  primary  arrivals  at  the  depth  of  the  layer.  A  near  vertical  reflection  for  a  deeper 
layer  would  have  an  arrival  shape  much  like  that  of  a  sideways  travelling  (backscattered) 
body  wave,  but  there  are  no  reflectors  below  the  seafloor  in  this  model.  S-wave  arrivals 
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Figure  3.8.  Arrivals  expected  at  a  vertical  array  of  geophones. 
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Figure  3.9.  Vertical  seismic  profile  (VSP)  for  model  SIN360.  Geophone  receivers  are 
located  at  zero  range  and  down  to  1.48  km  with  a  depth  spacing  of  40  meters.  Vertical 
displacement  is  shown.  Backscattered  shear  waves  are  marked. 
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would  be  similar  but  with  a  different  particle  motion  than  P-wave  arrivals. 

The  seismograms  of  vertical  particle  motion  are  shown  in  figure  3.9  for  the  vertical 
array  of  geophones  at  zero  range.  At  shallow  depths,  the  high  amplitude  Stoneley  waves 
dominate  the  signal  after  the  direct  wave.  However,  the  interface  waves  only  have  an  effect 
on  the  first  few  receiver  locations  so  that  at  depth,  only  body  waves  appear  on  the 
seismograms.  Because  of  the  fact  that  the  interface  waves  die  off  so  rapidly,  the  vertical 
geometry  is  very  well  suited  for  looking  at  both  primary  and  scattered  body  waves.  The 
primary  compressional  wave  is  the  first  arrival  on  both  plots.  Primary  shear  waves  are  not 
converted  at  the  interface  for  normal  incidence  compressional  energy.  Since  there  are  no 
sharp  breaks  in  the  velocity  profile  of  the  model  (other  than  the  seafloor)  all  of  the  rest  of 
the  energy  in  the  seismograms  must  be  scattered  body  waves. 

It  can  be  assumed  from  the  results  seen  in  the  snapshots  that  most  of  the  scattered 
body  wave  energy  is  in  the  form  of  shear  waves.  Since  a  large  part  of  the  energy  is 
travelling  nearly  horizontally  (at  least  at  later  times)  the  relatively  large  amplitudes  in  the 
vertical  motion  seismogram  would  also  indicate  scattered  shear  waves.  The  shapes  of  the 
scattered  arrivals,  if  they  can  be  singled  out,  are  indicative  of  the  curvature  of  the  wavefront 
and  the  direction  (relative  to  horizontal)  it  is  travelling.  Simplistically,  if  a  scattered  plane 
wave  is  travelling  upwards,  it  will  have  a  negative  slope  in  the  t-x  space  of  figure  3.9.  A 
downward  travelling  wave  will  have  positive  slope  in  the  seismogram.  Curvature  of  the 
wavefront  will  cause  a  slight  curvature  in  the  seismogram.  Arrivals  in  figure  3.9  which  are 
curved  and  have  a  generally  negative  slope  are  marked  as  backscattered  S-waves.  These 
are  not  the  only  backscattered  S-waves  but  are  pointed  out  here  because  this  type  of  arrival 
has  also  been  interpreted  as  being  caused  by  deep  reflectors  in  the  ocean  crust  (Becker,  et 
al„  1988;  Kennet  et.al.,  1980).  The  unique  shape  of  these  arrivals  on  a  vertical  array  is 
due  to  the  backscattered  energy  which  has  reached  depth  and  started  to  turn  back  towards 
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the  seafloor.  The  relatively  high  velocities  at  depth  causes  the  arrival  times  at  depth  to  be 
earlier  than  the  arrival  times  near  the  water-solid  interface.  S;nce  these  wavefronts  can  also 
be  clearly  seen  in  the  snapshots,  it  is  our  contention  that  these  arrivals  seen  in  vertical 
seismic  arrays  are  not  due  to  near  vertically  travelling  P-waves,  but  rather,  from  near 
horizontally  travelling  shear  waves  which  have  turned  at  depth. 

Three  other  initial  models  were  run  to  study  the  general  effects  of  the  influence  of 
sinusoid  wavelength  and  amplitude  on  the  scattering  of  seismic  energy.  The  seafloor 
profiles  of  these  models  are  shown  in  figure  3.2  and  have  wavelengths  of  900  meters  and 
1.8  kilometers.  Two  of  the  models,  SIN90A  and  SIN90B  have  the  same  wavelength  but 
different  amplitudes  (150  meters  and  80  meters  respectively).  The  different  seafloor 
topographic  profiles  are  manifest  in  two  ways.  First,  the  angle  of  incidence  vs.  range  plots 
for  the  three  models  (figure  3.10)  show  many  more  critical  angle  crossings  than  for 
SIN360  which  has  a  much  larger  wavelength.  Sharper  topography,  especially  in  model 
SIN90A  causes  much  more  complicated  horizontal  gradients  to  be  present  just  below  the 
seafloor. 

Results  for  model  SIN180  are  shown  in  figure  3.1 1.  Many  of  the  same  features 
seen  in  SEN360  are  also  present  in  this  model.  Additional  diving  waves  of  both 
compressional  and  shear  types  are  present  in  SIN  180  and  can  be  seen  in  both  the  4.0 
second  snapshot  and  time  series  plots.  Additional  P-waves  seen  in  the  time  series  were 
not  present  in  SIN360  and  are  due  to  the  slightly  more  variable  slopes  present  with  the 
shorter  wavelength  topography.  Also,  many  more  additional  S-waves  are  possible 
(seefigures  3.7a  and  3.7b)  because  of  the  steeper,  shorter  hills.  In  general,  more  energy  is 
able  to  enter  the  crust  at  large  ranges  as  shear  energy  than  in  SIN360  because  sections  of 
topography  with  angles  of  incidence  less  than  the  critical  angle  are  present  even  at  large 
ranges  (see  figures  3.10  and  3.1 1).  Of  course,  steeper,  shorter  topography  also  produces 
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Figure  3.10.  Angle  of  incidence  vs.  range  for  models  SIN90B  (3.10a.)  SIN90A.  (3.10b.) 
and  SIN180  (3.10c.).  Wavelength  of  the  topography  is  90  meters  for  SIN90A  and 
SIN90B  and  1.8  km.  for  SIN180.  Amplitude  for  the  topography  is  150  meters  for 
SIN90B  and  SIN  180  and  80  meters  for  SIN90A. 
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Figure  3.11.  Results  from  model  SIN  180.  3.11a.  The  4.0  second  snapshot  of 
compressional  (top)  and  shear  (bottom)  energy.  Note  the  presence  of  multiple  P  and  S 
diving  waves.  3. 1  lb,c.  Waterborne  receivers  (500  meters  above  the  seafloor  baseline)  for 
model  SIN  1 80  at  two  different  amplitude  scales.  3. 1  ld,e.  Pressure  time  series  for 
seafloor  receivers  in  model  SIN  180  at  two  different  amplitude  scales.  Multiple  P  and  S 
diving  waves  complicate  the  first  arrival  picks  for  these  phases. 
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Figure  3.11.  Results  from  model  SIN180.  3.11a.  The  4.0  second  snapshot  of 
compressional  (top)  and  shear  (bottom)  energy.  Note  the  presence  of  multiple  P  and  S 
diving  waves.  3.1  lb, c.  Waterborne  receivers  (500  meters  above  the  seafloor  baseline)  for 
model  SIN  180  at  two  different  amplitude  scales.  3.1  ld,e.  Pressure  time  series  for 
seafloor  receivers  in  model  SIN  180  at  two  different  amplitude  scales.  Multiple  P  and  S 
diving  waves  complicate  the  first  arrival  picks  for  these  phases. 
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Figure  3.11.  Results  from  model  SIN  180.  3.11a.  The  4.0  second  snapshot  of 
compressional  (top)  and  shear  (bottom)  energy.  Note  the  presence  of  multiple  P  and  S 
diving  waves.  3.1  lb,c.  Waterborne  receivers  (500  meters  above  the  seafloor  baseline)  for 
model  SIN  180  at  two  different  amplitude  scales.  3. 1  ld,e.  Pressure  time  series  for 
seafloor  receivers  in  model  SIN  180  at  two  different  amplitude  scales.  Multiple  P  and  S 
diving  waves  complicate  the  first  arrival  picks  for  these  phases. 
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more  dramatic  travel  time  anomalies  than  in  SIN360.  Also  note  that  the  scattered  noise  in 
the  seafloor  receivers  is  about  the  same  for  all  four  of  the  initial  models.  This  suggests  that 
it  is  probably  due  to  the  10  meters  steps  of  the  interface,  which  is  about  the  same  for  all  4 
models  (i.e.  they  have  approximately  the  same  micro-roughness). 

Models  SIN90A  and  SIN90B  were  chosen  not  only  to  demonstrate  shorter 
wavelength  features  but  also  to  compare  models  with  the  same  wavelength  but  different 
amplitudes.  The  snapshots  and  time  series  for  these  models  are  shown  in  figures  3.12  and 
3.13.  SIN90A  has  perhaps  the  most  obvious  effects  as  any  of  the  initial  models  because  of 
its  relatively  steep  facets  and  wide  fluctuations  in  incidence  angle  with  range.  Although  the 
angle  of  incidence  vs  range  curve  for  SIN90B  has  as  many  fluctuations,  they  are  not  as 
large  in  amplitude  as  in  SIN90A  (see  figure  3. 10).  The  principle  effect  of  the  fluctuations 
in  angle  of  incidence  is  the  presence  of  additional  P  and  S  diving  waves.  The  P-diving 
arrival  of  SIN90A  is  made  up  of  a  number  of  diving  waves  with  no  distinct  principle 
arrival.  The  shear  diving  wave  arrival  is  also  extremely  incoherent  and  practically 
indistinguishable  from  'noise'. 

The  compressional  and  shear  diving  waves  in  model  SIN90B  (figure  3.13)  are  both 
much  more  coherent  but  additional  waves  can  still  be  resolved.  The  shallower  topography 
of  SIN90B  causes  the  diving  wave  arrivals  to  appear  more  like  those  for  the  laterally 
homogeneous  FLAT  model  of  Dougherty  and  Stephen  (1987,1988;  chapters  1  and  2). 
Further  reduction  in  the  height  of  the  topography  would  cause  fewer  additional  diving 
waves  and  an  increase  in  the  coherence  of  the  primary  arrivals. 

Ray  theoretical  shadow  zones  occur  beyond  ranges  where  the  source  wave  is 
tangent  to  the  seafloor.  However,  with  the  full  wave  solution  of  the  finite  difference 
technique,  'shadow  zones'  are  impossible  because  of  diffracted  'Franz-type' 
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circumferential  waves  (Doolittle  et  al.,  1967;  Harbold  and  Steinberg,  1968;  Neubauer, 
1968).  Areas  in  ray  theory  shadow  zones  occur  at  ranges  greater  than  three  kilometers  in 
SIN90A  but  not  until  ranges  greater  than  nine  kilometers  in  SIN90B.  This  causes  the 
direct  and  reflected  arrivals  of  SIN90A  to  be  less  continuous  than  those  of  SIN90B. 
Backscattered  reflected  P-waves  will  occur  whenever  angle  of  incidence  is  equal  to  zero. 
These  account  for  the  first  few  strong  backscattered  arrivals  in  the  water  column  receivers 
of  the  two  models. 

The  seafloor  receivers  for  both  models  contain  bands  of  high  amplitude  interference 
patterns  above  valleys.  This  is  most  obvious  for  SIN90A  (figure  3.12b-e).  Since  these 
bands  do  not  show  up  in  the  water  column  receivers,  they  must  be  due  to  constructive 
interference  of  interface  waves  travelling  in  different  directions  over  the  topography.  A  less 
regular  seafloor  structure  would  not  tend  to  produce  these  regular  areas  of  interference 
because  of  less  coherent  interaction  of  the  wavefronts. 

A  number  of  important  scattering  phenomena  have  been  found  in  the  previous 
models.  It  appears,  however,  that  there  are  features  of  these  models  which  cannot  be  due 
to  scattering  from  the  larger  sinusoidal  topography  but  must  be  due  to  scattering  from 
smaller  scale  features  such  as  the  steps  comprising  the  sloping  seafloor.  Intuitively,  one 
can  identify  those  scattering  characteristics  which  are  due  to  larger  scale  sinusoidal  features 
(i.e.  multiple  diving  waves,  backscattered  body  waves  for  negative  angle  of  incidence,  and 
anomalies  in  travel  time  and  amplitude)  and  those  due  to  microroughness  of  the  sloping 
seafloor  (i.e.  interface  waves  and  the  majority  of  scattered  body  waves).  Also,  if 
scattering  from  ten  meter  steps  is  indeed  prominent,  then  it  probably  leads  more  to  Bragg 
scattering  problems  than  the  regular  nature  of  the  sine  wave  (at  least  for  the  larger 
wavelength  models).  Fortunately,  the  superposition  of  ten  meter  features  on  larger  scale 
topography  is  probably  not  unrealistic  for  many  parts  of  the  seafloor,  particularly  those 
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Figure  3.12.  Results  from  model  SIN90 A.  3.12a.  The  4.0  second  snapshot  of 
compressional  (top)  and  shear  (bottom)  energy.  Note  the  presence  of  multiple  P  and  S 
diving  waves.  3.12b,c.  Waterborne  receivers  (500  meters  above  the  seafloor  baseline)  for 
model  SIN90A  at  two  different  amplitude  scales.  3.12d,e.  Pressure  time  series  for 
seafloor  receivers  in  model  SIN90A  at  two  different  amplitude  scales. 
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Figure  3. 12.  Results  from  model  SIN90A.  3.12a.  The  4.0  second  snapshot  of 
compressional  (top)  and  shear  (bottom)  energy.  Note  the  presence  of  multiple  P  and  S 
diving  waves.  3.12b,c.  Waterborne  receivers  (500  meters  above  the  seafloor  baseline)  for 
model  SIN90A  at  two  different  amplitude  scales.  3.12d,e.  Pressure  time  series  for 
seafloor  receivers  in  model  SIN90A  at  two  different  amplitude  scales. 
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Figure  3.12.  Results  from  model  SIN90A.  3.i2a.  The  4.0  second  snapshot  of 
compressional  (top)  and  shear  (bottom)  energy.  Note  the  presence  of  multiple  P  and  S 
diving  waves.  3.12b,c.  Waterborne  receivers  (500  meters  above  the  seafloor  baseline)  for 
model  SIN90A  at  two  different  amplitude  scales.  3.12d,e.  Pressure  time  series  for 
seafloor  receivers  in  model  SIN90A  at  two  different  amplitude  scales. 
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Figure  3.13.  Results  from  model  SIN90B.  3.13a.  The  4.0  second  snapshot  of 
compressional  (top)  and  shear  (bottom)  energy.  Note  the  presence  of  multiple  P  and  S 
diving  waves.  3.13b,c.  Waterborne  receivers  (500  meters  above  the  seafloor  baseline)  for 
model  SIN90B  at  two  different  amplitude  scales.  3.13d,e.  Pressure  time  series  for 
seafloor  receivers  in  model  S 
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Figure  3.13.  Results  from  model  SIN90B.  3.13a.  The  4.0  second  snapshot  of 
compressional  (top)  and  shear  (bottom)  energy.  Note  the  presence  of  multiple  P  and  S 
diving  waves.  3.13b,c.  Waterborne  receivers  (500  meters  above  the  seafloor  baseline)  for 
model  SIN90B  at  two  different  amplitude  scales.  3.13d,e.  Pressure  time  series  for 
seafloor  receivers  in  model  S 
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Figure  3.13.  Results  from  model  SIN90B.  3.13a.  The  4.0  second  snapshot  of 
compressional  (top)  and  shear  (bottom)  energy.  Note  the  presence  of  multiple  P  and  S 
diving  waves.  3.13b,c.  Waterborne  receivers  (500  meters  above  the  seafloor  baseline)  for 
model  S1N90B  at  two  different  amplitude  scales.  3.13d,e.  Pressure  time  series  for 
seafloor  receivers  in  model  S 
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areas  close  to  mid  ocean  ridges.  Therefore,  scattering  phenomena  due  to  this 
microroughness  are  not  unexpected  in  these  areas. 

The  characteristic  criss-cross'  pattern  of  scattering  from  small  scale  roughness  on 
interfaces  in  this  work  is  similar  to  that  seen  by  Paul  and  Campillo  (1988)  in  their  work  on 
the  diffraction  of  elastic  waves  from  corrugated  surfaces.  It  is  unclear  from  their  results 
whether  the  diffracted  energy  is  due  entirely  to  the  'real'  corrugations  or  if  some  of  it  is 
indeed  due  to  the  discretization  of  the  interface.  Hill  and  Levander  (1984 )  and  Levander 
and  Hill  (1985)  also  observed  this  pattern  of  scattering  in  finite  difference  models  with 
small  scale  rectangular  roughness.  Levander  (1988)  also  presents  the  results  from  a  fourth- 
order  finite  difference  scheme  with  a  smoothly  varying  seafloor  model  (seafloor 
topography  characterized  by  one  cycle  of  a  cosinusoid).  In  this  model  there  is  a  large 
amount  of  scattered  energy  which  is  attributed  solely  to  scattering  from  the  large  scale 
sinusoid  feature.  However,  this  feature  is  probably  too  large  (topography  amplitudes  00 
m.,  wavelengths  1.0  km,  and  seismic  frequency=30  Hz)  with  respect  to  the  seismic 
wavelength  to  account  for  all  of  this  energy.  Some,  if  not  most  of  this  scattered  energy 
must  be  due  to  the  relatively  large  (10  meter)  grid  steps  used  to  define  the  interface. 
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TEST  MODELS 


The  remainder  of  this  chapter  will  be  devoted  to  the  determination  of  the  validity 
and  accuracy  of  scattering  from  the  microroughness  caused  by  the  grid  discretization  of 
rough  interfaces.  Two  tests  were  devised  to  investigate  the  accuracy  of  the  finite  difference 
solution  to  the  rough  seafloor  problems.  First,  the  grid  increments  were  decreased  in  order 
to  increase  accuracy  by  allowing  more  grid  points  per  wavelength  and  to  more  accurately 
describe  a  smoothly  varying  interface.  The  standard  15  points  per  wavelength  used  for  the 
preliminary  models  is  probably  not  fine  enough  to  solve  the  sinusoidal  model  but  may  be 
sufficient  to  solve  the  actual  stepwise  interface  imposed  by  the  grid.  The  second  way  to 
check  accuracy  of  the  solution  is  by  solving  the  reciprocal  problem.  A  solution  to  the 
elastodynamic  system  should  be  reciprocal  in  both  space  and  time  (Aki  and  Richards,  1980 
p.25-28).  Two  r ...  ..s  with  geometries  which  are  reciprocal  in  space  are  presented. 

Fine'*  grids 

The  finite  difference  technique  is  a  very  computationally  intensive  method  for  the 
solution  of  the  elastic  wave  equation.  At  the  standard  stability  condition  of  ten  points  per 
wavelength,  the  computer  power  exists  to  solve  realistic  seismic  refraction  models 
relatively  easily.  Changing  grid  parameters  can  quickly  increase  the  size  of  the  model 
beyond  the  limits  of  an  economic  solution.  Halving  the  grid  step  size  (from  ten  to  five 
meters)  effectively  increases  the  size  of  the  model  eight  times.  The  number  of  grid  steps  in 
range  and  depth  must  both  be  doubled  to  result  in  a  model  with  the  same  real  range  and 
depth  in  kilometers.  In  addition,  the  time  step  must  also  be  halved  (and  the  number  of 
steps  doubled)  to  preserve  the  stability  condition  in  time  (Dougherty  and  Stephen,  1988 
(chapter  2);  Kelly  etai.,  1976;  Virieux,  1986).  The  computer  resources  available  during 
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this  work  allowed  the  grid  dimension  to  be  halved  to  five  meters  per  step  but  did  not  allow 
for  any  grids  at  finer  increments. 

Two  models  were  computed  for  the  finer  grid  increment  of  five  meters.  Both  of 
these  models  had  the  same  large  scale  sinusoidal  topography  and  velocity  profiles  as  model 
SIN360.  The  first  of  these  models,  TEST1,  used  the  same  interpolation  scheme  to  define 
the  seafloor  as  in  SIN360  but  the  smaller  grid  steps  allowed  a  smoother  interface  to  be 
defined  (see  figure  3. 14).  The  small  five  meter  steps  should  have  a  lesser  effect  on  the  10 
Hz  seismic  waves  than  did  the  ten  meter  steps  of  SIN360.  The  second  model,  TEST2,  has 
the  same  ten  meter  step  structure  of  the  water-solid  interface  of  SIN360  (see  figure  3.14). 
However,  the  five  meter  spacing  should  provide  a  more  accurate  solution  than  the  coarser 
ten  meter  grid. 

TEST1  was  slightly  shorter  in  range  than  the  initial  models  and  has  overall  grid 
dimensions  of  1900  by  550  grid  points  (9.5  x  2.75  kilometers  at  five  meters  per  grid 
point)  with  the  same  boundary  and  source  conditions  of  the  initial  models  (see  figure  3.1). 
The  distance  between  the  bottom  of  the  transition  zone  and  the  top  of  the  telegraph  region  is 
slightly  less  for  model  TEST1  but  this  should  not  have  a  significant  effect  since  we  are 
interested  mainly  in  interface  effects.  The  model  was  run  out  to  8100  time  steps  (4.05 
seconds  at  0.0005  seconds/time  step)  to  facilitate  comparison  of  the  4.0  second  snapshot  of 
compressional  and  shear  energy  between  the  test  models  and  SIN360.  This  snapshot  for 
TEST1  is  shown  in  figure  3.15a.  In  general,  the  amount  of  scattered  energy  is  much  less 
in  TEST1  when  compared  to  the  4.0  second  snapshot  of  SIN360  (figure  3.60-  A  few 
isolated  scattered  shear  waves  are  still  present  in  TEST1  but  are  of  a  much  lesser  energy 
level  than  most  of  the  corresponding  waves  in  SIN360.  Also,  the  Stoneley  waves  at  the 
interface  are  much  lower  in  amplitude  than  in  the  previous  models.  Some  scattering 
phenomena  which  are  presumably  due  to  the  effect  of  the  larger  scale  sinusoidal 
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Figure  3.14.  Close-up  of  the  finite  difference  grid  interfaces  for  models  SIN360  (top), 
TEST1  (middle),  and  TEST2  (bottom).  The  grid  interface  is  shown  by  the  solid  line  with 
x's  (x's  being  located  at  grid  points).  The  actual  sinusoidal  boundary  is  shown  by  the  plain 
solid  line.  Grid  spacing  is  10  meters  for  SIN360  and  5  meters  for  models  TEST1  and 
TEST2. 
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Figure  3.15.  Results  from  model  TEST  1.  This  model  has  the  same  large  scale  seafloor 
topography  as  model  SIN360  (figures  3.6  and  3.7)  but  with  a  finer  grid  discretization. 

This  results  in  a  smoother  definition  of  the  seafloor  on  the  finite  difference  grid  (see  figure 
3.14).  3.15a.  The  4.0  second  snapshot  of  compressional  (top)  and  shear  (bottom)  energy. 
Note  that  much,  but  not  all,  of  the  incoherent  scattered  'noise'  has  not  been  generated  in 
this  model.  Multiple  S-diving  waves  however,  are  still  present. 
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rigure  3.15b.  Results  from  model  TEST1.  This  model  has  the  same  large  scale  seafloor 
topography  as  model  SIN360  (figures  3.6  and  3.7)  but  with  a  finer  grid  discretization. 
This  results  in  a  smoother  definition  of  the  seafloor  on  the  finite  difference  grid  (see  figure 
3.14).  Presented  here  are  the  time  series  for  seafloor  receivers.  Much,  but  not  all  of  the 
incoherent  scattered  energy  following  the  direct  and  P-diving  waves  has  not  been 
generated.  Some  low  amplitude  backscattered  interface  wave  energy  is  still  present 
following  the  direct  wave. 
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topography  are  still  present  in  TEST1  as  expected.  These  include  additional  shear  diving 
waves,  and  travel  time  and  amplitude  anomalies  of  the  principle  arrivals  (when  compared  to 
the  FLAT  model). 

The  same  general  results  are  also  evident  in  the  seafloor  time  series  of  TEST1 
(figure  3.15b).  The  amplitude  of  the  'noise'  following  the  direct  wave  is  much  reduced  but 
there  are  still  a  few,  low  amplitude,  scattered  shear  wave  arrivals  present  in  the  time  series. 
Interface  wave  amplitude  following  the  direct  wave  is  less  with  the  finer  step  structure  of 
TEST1  as  well.  Also,  as  in  the  snapshots,  characteristics  due  to  the  larger  scale 
topography  are  present  and  show  up  more  clearly  than  in  SIN360.  Since  only  seafloor 
receiver  time  series  were  saved  for  this  model,  the  signal  from  higher  in  the  water  column 
is  not  presented. 

The  results  of  TEST  1  indicate  that  the  majority  of  the  backscattered  energy  in 
SIN360  cannot  be  due  to  interactions  with  the  large  scale  topography.  If  this  scattering  is 
due  to  diffractions  from  the  smaller  ten  meter  steps  on  the  finite  difference  seafloor,  then 
what  remains  to  be  determined  is  if  this  is  a  realistic  phenomenon  at  10  Hz  wavelengths,  or 
is  it  due  to  inaccuracies  of  the  coarse  grid.  The  ideal  procedure  for  obtaining  an  answer  to 
this  problem  would  be  to  calculate  models  with  decreasing  grid  size  until  a  constant 
solution  is  obtained.  As  mentioned  earlier,  the  computer  resources  available  limited  this 
procedure  to  only  one  iteration  with  smaller  grid  increments  (five  meter  spacing  instead  of 
10  meter).  However,  useful  conclusions  can  be  drawn  even  from  this  limited  test. 

TEST2  used  an  interpolation  scheme  for  seafloor  definition  which  created  ten  meter 
steps  on  the  five  meter  grid  which  were  identical  to  the  ten  meter  steps  of  SIN360  (see 
figure  3.14).  Using  this  scheme,  the  diffractions  from  the  coarse  seafloor  can  be  modelled 
more  accurately  on  the  finer  grid.  TEST2  has  the  same  depth  dimension  as  TEST1  (550 
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points)  but  only  reaches  a  range  of  6.0  kilometers  (  1200  points  at  five  meters  per  point). 
The  model  size  is  more  than  sufficient  for  close  range  ( less  than  five  kilometers)  scattering. 

The  4.0  second  snapshot  of  compressional  and  shear  energy  and  the  seafloor 
pressure  time  series  for  model  TEST2  are  shown  in  figures  3.16a  and  3.16b,  respectively. 
The  best  qualitative  comparison  can  be  made  with  the  4.0  second  snapshots  of  TEST2 
(figure  3.16a)  and  SIN360  (figure  3.6f).  In  general,  the  strength  of  scattered  energy  is 
about  the  same  in  the  two  models.  In  both  cases,  scattering  from  the  stepwise  seafloor  is 
much  stronger  that  in  the  TEST1  model  which  has  a  smoother  grid  defined  seafloor.  A 
comparison  between  the  seafloor  pressure  time  series  of  models  SIN360  (figure  3.7)  and 
TEST2  (figure  3.16b.)  also  shows  general  agreement  between  scattered  amplitudes. 
However,  because  TEST2  represents  a  more  accurate  solution  to  the  problem  (because  of 
finer  grid  spacing),  corresponding  individual  traces  from  the  two  models  are  not  identical. 
Results  from  both  the  snapshots  and  time  series  indicate  that  features  on  the  order  of  1/15 
wavelength  produce  an  important  component  of  the  scattered  field. 

Reciprocal  problem 

A  second  way  to  verify  the  accuracy  of  the  finite  difference  formulation  used  is  to 
compute  the  reciprocal  problem.  Because  elastic  wave  propagation  is  a  linear  process,  the 
signal  generated  between  a  source-receiver  pair  should  be  independent  of  direction.  That 
is,  the  process  of  elastic  wave  propagtion  is  reciprocal  in  the  spatial  dimension  (Aki  and 
Richards,  1980,  p.25-28).  For  example,  within  the  same  environment  or  model,  the 
pressure  signal  received  by  a  receiver  at  a  range  of  2.0  kilometers  from  a  source  at  0.0 
kilometers  should  be  the  same  as  a  signal  received  by  a  pressure  receiver  at  0.0  kilometers 
range  from  a  source  at  2.0  kilometers  (for  source  and  receivers  at  reciprocal  heights  and  in 
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Figure  3.16.  Results  from  model  TEST2.  This  model  has  the  same  large  scale  sinusoidal 
seafloor  structure  as  models  SIN360  and  TEST1.  However,  the  finer  discretization  of  the 
seafloor  in  model  TEST2  follows  the  coarse  steps  of  model  SIN360  (see  figure  3. 14). 
3.16a.  The  4.0  second  snapshot  of  compressional  (top)  and  shear  (bottom)  energy. 
Incoherent  scattered  energy  (scattering  from  stepwise  seafloor  discretization)  is  generated 
and  is  as  strong  as  in  model  SIN360  (see  figure  3.6). 
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Figure  3. 16b.  Results  from  model  TEST2.  This  model  has  the  same  large  scale  sinusoidal 
seafloor  structure  as  models  SIN360  and  TEST1.  However,  the  finer  discretization  of  the 
seafloor  in  model  TEST2  follows  the  coarse  steps  of  model  SIN360  (see  figure  3.14). 
Presented  here  are  the  pressure  time  series  for  seafloor  receivers.  Scattering  from  the 
stepwise  seafloor  is  as  prevalent  as  in  model  SIN360  (see  figure  3.7). 
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a  fluid). 


In  order  to  calculate  this  problem  using  our  formulation,  a  few  changes  had  to  be 
made  on  the  grid  used  for  all  of  the  previous  models.  Figure  3. 17  shows  the  grid 
geometry  used  for  these  two  reciprocal  models.  It  is  important  that  the  actual  source  and 
the  receivers  be  at  reciprocal  heights  in  the  water  column.  In  this  case  they  are  at  the  same 
height  because  this  formulation  does  not  allow  for  placing  the  source  too  near  the  water- 
solid  interface.  The  source  and  receivers  are  located  1.5  kilometers  above  the  topography 
baseline.  This  was  used  as  the  optimal  level  location  to  minimize  false  reflections  from  the 
top  and  comers  of  the  finite  difference  grid.  Also,  to  help  minimize  false  reflections,  the 
entire  water  column  (three  kilometers  thickness)  must  be  calculated  for  these  models  which 
dramatically  increases  computation  time.  Only  the  bottom  kilometer  of  the  water  column 
was  calculated  using  finite  differences  in  the  previous  models  with  the  rest  of  the  energy 
being  absorbed  by  the  top  boundary.  This  top  absorbing  boundary  was  also  used  for  the 
reciprocal  problem  but  since  it  is  not  perfectly  absorbing,  it  was  necessary  to  place  it  as  far 
away  from  the  receivers  as  possible.  Since  the  source  is  actually  introduced  as  downward 
travelling  energy  along  a  line  near  the  seafloor  (see  figure  3. 13)  its  direct  signal  does  not 
appear  in  the  time  series  (of  receivers  higher  up  in  the  water  column). 

Two  models  were  run  (RECIP1  and  RECIP2)  which  had  reciprocal  structure  and 
geometry.  A  grid  dimension  of  10  meters  was  used  for  both  models  with  a  time  step  of 
0.001  second.  A  constant  velocity  crust  (3.2  km/sec  for  P-waves,  1.85  km/sec  for  S- 
waves)  beneath  the  sinusoidal  topography  was  used  rather  than  the  gradients  and  lateral 
heterogeneity  profiles  of  the  previous  models.  This  was  done  to  isolate  effects  of  the  small 
interface  steps  from  any  energy  turning  in  the  gradients.  Both  models  used  the  same 
interface  structure  as  SIN360  (topography  amplitude  =  150  meters,  wavelength=3.6 
kilometers)  with  receivers  located  every  40  meters.  Although  only  a  few  receivers  will  be 
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used  for  the  reciprocity  results,  the  entire  line  was  saved  and  used  for  identification  of 
arrivals.  The  topography  of  RECIP2  used  the  reciprocal  sinusoid  to  RECIP1  (i.e.  a  phase 
shift  of  180  degrees  was  added,  see  figure  3.17),  The  periodic  nature  of  the  seafloor  in 
both  models  allows  more  than  one  source-receiver  pair  in  each  model  to  overlap  with  the 
corresponding  pair  in  the  other  model.  RECIP1,  with  its  source  over  a  valley  and  receivers 
at  1.8  and  5.4  kilometers  (and  so  on  every  3.6  km.  in  range)  over  hills  is  reciprocal  to 
RECIP2,  with  its  source  over  a  hill  and  receivers  at  1.8  and  5.4  kilometers  (and  out)  over 
valleys.  The  signals  from  the  1.8  kilometer  receivers  of  the  two  models  should  be  the 
same,  as  well  as  those  for  the  5.4  kilometer  receivers. 

It  is  important  to  understand  that  the  time  series  from  all  of  the  corresponding 
receivers  in  the  two  models  will  not  be  similar,  only  those  source-receiver  pairs  which  have 
reciprocal  geometry.  Figure  3.18  presents  the  time  series  signals  from  two  reciprocal 
receiver  pairs  at  1.8  and  5.4  kilometers  in  models  RECIP1  and  RECIP2.  In  both  cases, 
the  agreement  is  very  good,  especially  at  earlier  times  in  the  time  series.  Slight 
discrepancies  at  later  times  are  due  to  two  factors.  First,  it  has  already  been  shown  that  the 
10  meter  solution  is  slightly  inaccurate  when  compared  to  the  5  meter  solution  to  the  rough 
interface  problem.  Therefore,  some  directional  inaccuracies  may  be  present  in  the 
reciprocal  models.  A  larger  source  of  disagreement  in  the  later  sections  of  the  time  series  is 
due  to  boundary  problems.  Although  the  absorbing  boundaries  are  considered  to  be  quite 
efficient  in  absorbing  energy,  they  are  not  perfect  and  add  to  the  signal  at  late  times.  The 
receivers  were  placed  as  far  away  as  possible  from  any  problem  areas  but  still  received 
some  boundary  noise.  Since  the  geometry  which  includes  the  seafloor,  source-receiver 
pair,  and  grid  boundaries  is  not  reciprocal,  any  comer  or  edge  reflections  will  not 
correspond  in  the  two  models. 
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Figure  3. 17.  Grid  layout  for  the  reciprocal  problem.  The  entire  water  column  is  calculated 
using  finite  differences  for  this  problem.  Crust  below  the  seafloor  is  of  a  constant  velocity 
to  isolate  seafloor  scattering  effects.  Note  that  while  the  source- seafloor-receiver 
geometries  in  the  two  models  are  reciprocal,  the  relationship  between  the  source-receivers 
and  the  artificial  numerical  boundaries  and  comers  is  not.  Reciprocal  traces  occur  at  1.8 
and  5.4  km.  range  (  and  so  on  every  3.6  km.). 
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Figure  3.18.  Reciprocal  traces  at  1.8  and  5.4  km.  for  the  two  reciprocal  models.  The 
initial  part  of  the  traces  at  both  ranges  match  exactly.  Most  of  the  discrepancy  at  later  times 
(especially  for  the  1 .8  km.  traces)  is  due  to  boundary  interaction  which  is  not  reciprocal. 
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The  results  from  the  test  models  point  out  an  important  numerical  result  of  this 
work.  That  is,  the  need  to  consider  more  than  numerical  stability  and  dispersion  when 
determining  the  necessary  grid  spacing.  Since  features  with  sizes  on  the  order  of  1/15  of 
the  seismic  wavelength  have  been  shown  to  cause  significant  scattering,  the  ability  of  the 
grid  to  accurately  represent  the  model  must  also  be  considered.  While  the  models  presented 
in  this  work  are  numerically  stable,  they  do  not  accurately  represent  the  problem  of 
scattering  from  the  larger  scale  sinusoidal  topography. 

Other  numerical  methods  such  as  the  pseudo-spectral  method  (Fomberg,  1987; 
Gazdag,  1981;  Kosloff  et  al.,  1984)  and  higher  order  finite  differences  (Bayliss  et  al., 

1986;  Dablain,  1986;  Levander,  1988;  Vidale  et  al.,  1985)  can  be  more  efficient  for  certain 
problems  than  the  second  order  finite  difference  scheme  used  for  this  work  because  fewer 
grid  points  per  wavelength  are  necessary  for  numerical  stability.  However,  these  methods 
will  also  be  sensitive  to  the  problem  of  the  definition  of  non-planar,  non-level  topography 
with  a  rectangular  grid.  Therefore,  as  with  second  order  formulations,  these  other  methods 
will  need  much  denser  grids  than  needed  for  numerical  stability.  (Fomberg,  1988)  has 
shown  that  this  problem  can  be  avoided  with  the  pseudo-spectral  method  for  simple  rough 
interfaces  by  deforming  the  grid  around  the  interface  so  that  the  interface  lies  directly  on  a 
grid  line  (no  'stair-stepping’  needed  for  non-level  interfaces).  Some  of  the  efficiency  of 
that  method  is  lost  by  using  deformed  grids  but  the  method  shows  good  results  for  simple 
rough  interfaces  between  homogeneous  media.  However,  logical  deformation  of  grids  for 
more  complicated  2-D  structures  and  velocity  profiles  is  probably  not  realistic. 
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CONCLUSIONS 


The  finite  difference  models  of  scattering  from  sinusoidal  seafloor  topography 
presented  in  this  paper  have  demonstrated  the  effects  of  scattering  from  both  macro-  and 
micro-roughness,  as  well  as  some  important  constraints  on  the  method  used.  The  initial 
models,  run  with  relatively  coarse  grids  show  that  effects  such  as  travel  time  anomalies, 
additional  compressional  and  shear  diving  waves  and  some  strong  back  reflected  arrivals 
are  due  to  the  rough  topography.  Steep  topography  allows  energy,  especially  shear 
energy,  to  enter  the  seafloor  even  at  great  ranges.  Ray  theoretical  shadow  zones  are  not 
present  because  of  Franz-type  waves  diffracted  into  areas  where  the  grazing  angle  is  less 
than  zero. 

It  was  suspected  from  the  initial  model  results  that  not  all  of  the  scattered  signal 
could  be  due  to  the  large  scale  sinusoidal  topography.  Phenomena  such  as  backscattered 
body  waves  (predominandy  shear  waves),  scattering  into  interface  waves,  and  strong 
interference  effects  are  most  likely  due  to  the  grid  imposed  micro-roughness  superimposed 
onto  the  larger  scale  topography.  Seen  in  a  vertical  seismic  profile  geometry,  the 
backscattered,  near  horizontally  travelling  S-waves  appear  very  similar  to  arrivals  which 
have  previously  been  interpreted  as  near  vertically  travelling  P-waves  reflected  from  layers 
deep  m  the  crust  All  of  these  scattered  wave  types  are  important  because  of  the  abundance 
of  smaller  scale  features  on  the  sea  floor  especially  near  mid-ocean  ridges. 

Models  with  finer  'jrids  and  both  smooth  and  coarse  water-solid  interfaces 
demonstrated  that  scattering  from  interface  steps  on  the  order  of  1/15  wavelength  is  realistic 
and  important.  However,  even  finer  grids  (>30  points/wavelength)  are  necessary  for 
better  definition  of  a  smooth  seafloor  (without  significant  micro-roughness).  The 
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reciprocal  problem  computed  with  a  coarse  grid  (15  points/wavelength)  demonstrated  that 
the  coarse  grid  provides  a  realistic  simulation  of  wave  propagation  through  sinusoidal 
topography  with  grid  imposed  small  scale  topography. 
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CHAPTER  4 


A  time  domain  solution  to  acoustic  wave  scattering  from  an 

infinite  elastic  cylinder 
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A  time  domain  solution  to  acoustic  wave  scattering  from  an 

infinite  elastic  cylinder. 


INTRODUCTION 


This  chapter  presents  a  time  domain  solution  to  the  problem  of  acoustic  wave 
scattering  from  an  infinite  elastic  cylinder  by  using  the  finite  difference  method.  It  is 
difficult  to  imagine  how  such  a  conceptually  'simple'  problem  has  been  the  subject  of  so 
much  interest  and  elucidation  over  the  years.  This  is  a  very  well  studied  problem  and  has 
been  so  for  over  100  years,  beginning  with  the  efforts  of  the  eminent  nineteenth  century 
acoustician,  Lord  Rayleigh  (Rayleigh,  1878).  Our  purpose  for  trying  this  problem  is  not  to 
shed  any  new  light  onto  the  mechanisms  of  scattering  or  to  discover  any  heretofore 
unknown  scattered  modes,  but  rather,  to  use  the  problem  as  a  stringent  test  of  the  finite 
difference  method  used  for  other,  more  geologically  relevant  scattering  problems 
(Dougherty  and  Stephen,  1988  (chapter  2);  Stephen,  1988;  Virieux,  1984;  Virieux,  1986). 

In  a  time  domain  sense,  two  phenomena  characterize  the  problem  of  scattering  from 
a  cylinder.  First,  the  backscattered  signal  received  after  the  illumination  of  an  infinite 
cylinder  by  an  incident  pulse  source  contains  not  a  single  reflected  arrival  but  a  number  of 
pulses  arriving  after  the  initial  specular  reflection.  Also,  insonification  by  a  CW 
(continuous  wavelet)  source  results  in  a  complex  interference  pattern  of  energy  scattered  in 
a  number  of  discrete  directions  from  the  cylinder.  Lord  Rayleigh  first  studied  this  problem 
for  the  cases  of  rigid  and  acoustic  cylinders  in  a  fluid  medium  and  found  that  much  of  the 
reverberated  energy  was  due  to  the  existence  of  diffracted  waves  which  travel  along  the 
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circumference  of  the  cylinder.  These  results  were  expanded  to  the  elastic  cylinder  case  by 
Far  an  (1951)  and  Lowan  (1946)  who  found  that  even  more  circumferential  wave  modes 
are  possible  when  the  cylinder  also  supports  shear.  The  situation  for  nonrigid  and  elastic 
cylinders  is  complicated  further  by  the  presence  of  transmitted  and  whispering  gallery 
modes  as  seen  by  Doolittle  et  al.  (1966,1968)  and  Brill  and  Uberall  (1970). 

Solutions  to  the  cylindrical  scattering  problem  have  been  obtained  by  a  number  of 
analytical,  experimental,  and  numerical  methods.  Analytically,  the  problem  has  been 
solved  by  two  basic  techniques  for  cylinders  with  radii  near  the  size  of  the  acoustic 
wavelength  or  larger.  Rayleigh  (1878),  Lowan  et  al.  (1946),  and  Faran  (1951)  presented 
solutions  which  were  comprised  of  series  of  normal  modes  of  the  cylinder.  This  method 
produces  a  solution  which  implicitly  contains  all  possible  wavetypes  scattered  by  the 
cylinder  if  enough  modes  are  included  The  disadvantage  of  this  method  is  that  individual 
wave  types  which  contribute  to  the  total  scattered  field  cannot  easily  be  individually 
examined.  Also,  a  fair  amount  of  computation  is  involved  in  making  the  series  converge. 
Other  authors,  (Doollittle  et  al.,  1966,1967,  Brill  and  Uberall,  1970,  and  Neubauer,  1986) 
have  employed  the  Sommerf eld- Watson  transformation  (Watson,  1919)  which  transforms 
the  infinite  normal  mode  series  expression  for  the  scattered  pressure  into  a  contour  integral 
in  the  complex  mode  number  plane.  While  this  solution  may  converge  more  quickly  than 
the  normal  mode  solution,  the  zeros  must  still  be  found  numerically.  The  disadvantages  of 
this  method  are  that  it  is  not  a  total  solution  (unless  all  of  the  zeros  are  found)  and  that  the 
numerical  search  for  the  poles  may  be  complicated,  especially  for  small  cylinder  radii  (Frisk 
et  al.,  1975;  Grace  and  Goodman,  1966).  Numerical  experiments  have  verified  theoretical 
results  (and  vice  versa)  for  both  the  total  field  solutions  and  for  individual  contributions  to 
the  scattered  field  (Brill  and  Uberall,  1970;  Faran,  1951;  Harbold  and  Steinberg,  1968; 
Neubauer,  1968;  Neubauer  et  al.,  1974;  Stoyanov  et  al.,  1989). 
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This  study  is  an  attempt  to  solve  the  cylindrical  scatterer  problem  by  using  the 
numerical  finite  difference  formulation  of  the  wave  equation.  Finite  difference  methods  are 
well  established  for  elastic  wave  propagation  problems  concerning  both  simple  and 
complex  structures  (Alterman  and  Loewenthal,  1972;  Dougherty  and  Stephen,  1988 
(chapter  2);  Frankel  and  Clayton,  1984;  Kelly  et  al.,  1976;  Levander  and  Hill,  1985; 
Stephen,  1988a,b;  Virieux,  1986).  The  robust  nature  of  the  method  is  not  always 
necessary  for  acoustics  problems  or  problem  with  simple  geometrical  relationships  (for 
which  analytical  solutions  are  easily  obtained).  While  various  finite  difference  formulations 
of  the  wave  equation  have  been  tested  against  analytical  solutions  for  layered  or  rectangular 
structures  (Kelly  et  al.,  1976;  Stephen,  1983;  Virieux,  1986),  the  more  rigorous  test  of 
scattering  from  non-planar,  non-rectangular  objects  has  not  been  well  studied.  The  finite 
difference  formulation  used  for  this  study  is  a  2-D  elastic  displacement-stress  code  in 
rectangular  coordinates  and  has  been  applied  to  earth  models  with  complex  structures 
(Dougherty  and  Stephen,  1988  (chapter  2);  Stephen,  1988;  Virieux,  1986).  The  models 
presented  in  this  work  are  solved  by  using  the  method  without  any  modifications  from  that 
used  for  the  earth  models.  Steady  state  (using  a  monochromatic  CW  source)  and  transient 
wave  (using  a  broad-band  Gaussian  pulse  source)  solutions  are  presented. 
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CW  SOURCE 


The  first  cylindrical  model  was  run  with  a  monochromatic  continuous  wavelet 
(CW)  source.  The  theory  for  this  situation  is  well  developed  for  a  homogeneous  elastic 
infinite  cylinder  and  is  briefly  outlined  below.  The  geometry  of  the  cylinder  and  incident 
CW  plane  wave  source  is  shown  in  figure  4.1  with  the  cylinder  extending  infinitely  in  both 
directions  along  the  Z  (cylindrical)  and/or  y  (Cartesian)  axes  and  the  plane  wave  extending 
infinitely  in  the  x  and  y  directions  (Cartesian).  The  normal  mode  solution  of  Faran  (1951) 
was  chosen  for  the  analytical  solution  because  of  its  relative  ease  of  calculation.  The 
Sommerfeld- Watson  transformation  technique  provides  insight  into  which  phases  are 
contributing  to  the  scattered  field  but  is  more  complicated  numerically  for  small  ka  values 
and  the  total  solution  is  not  easily  reached.  The  normal  mode  solution  gives  a  total 
scattered  pressure  v/avefield  which  is  easily  compared  to  the  total  field  produced  by  the 
finite  difference  method. 

The  solution  to  the  infinite  elastic  cylinder  problem  consists  of  solving  the  wave 
equation  inside  and  outside  of  the  cylinder  and  matching  solutions  at  the  boundary.  The 
infinite  cylinder  is  independent  of  the  Z-dimension  ( into  and  out  of  the  plane  of  the  page  in 
figure  4.2)  and  is  thus  able  to  be  modelled  easily  by  the  2-D  finite  difference  method  as 
well  as  by  the  analytical  methods.  Inside  the  homogeneous  elastic  cylinder,  displacements 
are  governed  by  the  elastic  wave  equation  for  homogeneous  media; 

r.2 

(1)  (  k  +  2^i)  V(  V*  u) -\i  Vx(  VxU)  =  p - 

5 1 
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Figure  4.1.  Geometrical  considerations  for  the  problem  of  plane  wave  scattering  by  an 
infinite  elastic  cylinder.  Cylindrical  coordinates  (bold)  were  used  for  theoretical 
formulations  and  Cartesian  coordinates  (in  parentheses)  were  used  for  the  finite  difference 
modelling.  The  cylinder  extends  infinitely  in  both  directions  along  the  Z  (cylindrical)  or  y 
(Cartesian)  axis. 
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Although  significant  ’coda'  exists  for  model  GAUSS  10,  suggesting  that  random 
scatter  is  important,  there  is  apparendy  little  change  in  the  waveform  caused  by  the  small 
scatterers.  Since  'ca  is  small  for  GAUSS  10,  the  wave  field  is  starting  to  see  the  velocity 
perturbations  as  an  equivalent  body  with  an  average  velocity  field  approaching  that  of  the 
FLAT  model.  The  random  noise  is  slightly  higher  in  frequencey  and  arises  from  scattering 
of  the  higher  frequency  components  of  the  source.  The  predominantly  10  Hz  P-diving 
wave  is  not  greatly  affected. 

The  distinction  between  the  cross  correlation  characteristics  of  the  different  models 
decreases  as  the  range  of  the  reference  trace  increases.  A  clear  separation  exists  between 
the  cross-correlation  patterns  of  the  different  models  if  the  reference  trace  is  at  a  relatively 
small  range  (such  as  6.0  km.  in  figure  2.17a.).  As  the  range  of  the  reference  trace 
increases  there  is  an  overall  decrease  in  most  of  the  absolute  values  of  the  correlation 
coeficients  as  well  as  a  decrease  in  the  distinction  of  the  patterns  between  models.  This  is 
probably  due  to  the  influence  of  the  interference  head  wave  at  greater  ranges.  That  is,  low 
correlation  between  traces  is  influenced  more  by  the  interference  head  wave  than  by  the 
heterogeneity  at  large  ranges.  Model  GAUSS50  which  has  a  lower  correlation  curve  than 
the  others  has  about  the  same  correlation  values  with  different  reference  ranges  because  the 
large  amount  of  scattering  causes  the  interference  head  wave  to  be  at  a  low  amplitude 
already  (see  figures  2.10  and  2.14  for  a  comparison  of  interference  head  wave  amplitudes). 
Models  with  less  random  scatter  have  stronger  interference  head  waves  at  large  ranges  so 
that  the  decrease  in  correlation  coefficients  with  range  is  larger  for  these  models. 
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Figure  4.2.  Grid  geometry  used  for  the  finite  difference  solution.  Absorbing  boundaries 
based  on  the  parabolic  equation  are  used  on  the  top  and  right  sides  of  the  grid.  The  left 
hand  side  is  an  axis  of  symmetry.  Energy  travelling  downward  is  attenuated  in  the  stipled 
region  by  using  the  telegraph  equation.  The  cylinder  is  placed  on  the  axis  of  symmetry  for 
the  CW  source  model  and  extends  infinitely  into  and  out  of  the  page. 
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The  field  in  the  fluid  outside  of  the  cylinder  is  somewhat  less  complicated  than  that 
within  the  elastic  medium.  Waves  within  the  fluid  obey  the  acoustic  wave  equation  for 
homogeneous  media; 


(7) 


V2p  =  ±li 
2  „  2 


c3  8  t 


where  p  is  pressure  and  C3  is  the  compressional  wave  velocity  in  the  fluid.  The  incident 
CW  plane  wave  can  be  represented  by; 

oo 

(8)  P,  =  PO  £  £n  (-j)"  Jn  (  ksT  )  COS  (  110  ) 

n=0 


where  po  is  amplitude  ,  £n  is  the  Neumann  factor  (£0=1;  £n=2,  n>0),  j=sqrt(-l),  and  k3  is 
the  wavenumber  of  acoustic  waves  in  the  fluid.  The  radial  component  of  the  CW  plane 
wave  source  is  all  that  is  needed  to  match  boundary  conditions,  it  is  given  by; 


(9)  Uu— Bt-ie.(-j)‘8(J“(%)C0,(ne)) 

P3  (0  n=0 

Because  of  symmetry  about  0=0,  the  outgoing  scattered  pressure  and  displacement  must  be 
given  by; 

oo 

(10)  Ps=£cn(Jn(k3r)-jNn(k3r  ))  COS  (  110  ) 

n=0 
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(11)  U 


1 


£  S((Jn(k3T)-jNn(k3r  ))cos(n0)) 

2-  c"  8  r 

p3CO'n=0  °r 

where  Nn  is  the  Bessel  function  of  the  second  kind  of  order  n  and  cn  are  the  coefficients 
needed  to  describe  the  scattered  pressure  field.  These  coefficients  are  found  by  matching 
boundary  conditions  across  the  cylinder  surface.  At  a  liquid-solid  interface  there  are  three 
conditions  which  must  be  met; 

i.  continuity  of  normal  stress  (in  the  fluid,  normal  stress  is  given  by  the  total 
pressure); 

(12)  Pi  ■*"  Ps  i 

ii.  continuity  of  the  normal  component  of  displacement; 

(13)  uiir+UStr=Ur 

iii.  tangential  stress  within  the  cylinder  must  vanish  at  the  interface  (liquid  will  not 
support  shear); 

(14)  Trt0=IrtZ=O 

Solving  these  simultaneous  equations  for  ps  is  complicated  because  of  the  infinite  series 
involved  and  the  number  of  constants  (cn)  which  must  be  determined.  This  has  been  done 
by  Faran  (1951)  who  arrived  at  an  expression  for  the  amplitude  of  scattered  pressure  in  the 
far  field  as  a  function  of  azimuth; 
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(15) 


-f-l 

k  k3  r 


£  £n  sin  Tin  exp  (  jTln  )  cos  (  n0  ) 

n=0 


The  expressions  for  the  scattering  phase  angles,  r\,  are  too  involved  for  presentation  here 
but  are  given  by  Faran  (1951)  and  Lowan  et  al.  (1946)  and  are  composed  of  complex 
expressions  of  Bessel  functions  and  their  derivatives. 

Equation  15  (including  the  expressions  for  Ti)  was  solved  out  to  15  terms  for 
comparison  with  the  finite  difference  runs  discussed  below.  Although  this  modal  solution 
is  analytical,  it  must  be  numerically  evaluated.  In  order  to  test  the  accuracy  of  the  numerical 
formulation  of  the  analytical  solution  used  for  this  work,  the  results  of  the  code  were 
compared  against  Faran's  results  for  metal  cylinders  with  excellent  agreement  Also, 
values  output  by  the  numerical  formulation  for  the  intermediate  scattering  angles,  rj,  were 
tested  against  tables  in  Faran  (1951)  and  Lowan  et  al.  (1946 ) ,  again  with  excellent 
agreement 

The  velocities  and  densities  of  metal  cylinders  presented  in  most  of  the  literature  are 
too  high  to  be  realistic  for  earth  materials  in  the  upper  crust.  The  metal  cylinders  commonly 
used  for  test  experiments  have  relatively  large  elastic  constants  (  a  p-wave  velocity  of  4.28 
km/sec,  shear  wave  velocity  of  2. 165  km/sec  and  a  density  of  8.5  gm/cc  for  brass).  These 
high  values  represent  a  velocity  contrast  of  over  200  per  cent  at  the  water-solid  interface 
and  are  cause  instabilities  in  the  finite  difference  modeling  scheme  used  here.  The  oceanic 
crustal  scattering  models  presented  in  Dougherty  and  Stephen  (1988)  and  heterogeneous 
crustal  halfspace  models  presented  by  Frankel  and  Clayton  (1984)  have  a  maximum  crustal 
velocity  variation  of  10  per  cent.  Since  this  study  is  an  attempt  to  verify  the  use  of  the  finite 
difference  formulation  for  scattering  in  the  upper  crust,  a  much  smaller  velocity  contrast 
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was  chosen  for  the  infinite  cylinders  used  in  this  work.  Because  most  effects  need  a  fairly 
significant  contrast  to  be  visible,  a  p-wave  velocity  of  2.6  km/sec,  s-wave  velocity  of  1 .5 
km/sec,  and  density  of  1.33  gm/cc  were  chosen  as  the  parameters  for  the  cylinder.  These 
elastic  parameters  provide  enough  contrast  to  show  effects  but  not  so  much  that  it  is 
unrealistic  for  the  earth  models  or  that  instabilities  in  the  formulation  will  be  developed. 

The  orientation  of  the  cylinder  on  the  finite  difference  grid  is  shown  in  figure  4.2. 
A  radius  of  160  meters  was  used  to  give  a  ka  value  (wavenumber,  k,  times  radius,  a)  of 
6.7  for  a  10  Hz  source.  An  attempt  was  made  to  use  the  same  finite  difference  formulation 
and  grid  setup  (including  the  same  boundary  conditions,  transition  zone  setup,  etc)  as  was 
used  for  the  earth  models  (Dougherty  and  Stephen,  1988  (chapter  2),  and  chapter  3). 

While  this  may  not  be  the  optimal  setup  for  solving  the  cylinder  problem,  we  did  not  feel 
that  it  was  necessary  or  valid  to  create  special  conditions  for  this  problem.  The  cylinder 
was  placed  on  the  axis  of  symmetry  in  order  to  save  computer  costs.  Since  this  is 
theoretically  a  symmetrical  problem  (with  a  symmetrical  solution),  it  is  not  necessary  to 
calculate  both  sides  of  the  cylinder.  Some  initial  models  were  run  which  calculated  the 
entire  cylinder  to  check  the  feasibility  of  putting  the  cylinder  on  the  axis  of  symmetry. 
These  results  indicated  an  inherent  asymmetry  in  the  way  the  grid  was  being  defined.  This 
problem  was  corrected  and  the  finite  difference  models  were  indeed  made  to  be  symmetric. 

Receivers  are  located  in  a  semi-circle  surrounding  the  center  of  the  cylinder  at  a 
distances  of  1.4  kilometers.  Azimuthal  distribution  of  receivers  is  one  per  degree.  A  10 
Hz  sinusoidal  CW  source  was  introduced  as  a  plane  wave  along  the  top  of  the  grid.  The 
source  has  only  a  vertical  component  of  displacement  so  that  the  solution  is  symmetric 
about  the  center  of  the  cylinder  (range=0).  Other  details  of  the  grid  such  as  source 
introduction,  absorbing  boundaries,  etc.  can  be  found  in  Dougherty  and  Stephen  (1988). 
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The  CW  finite  difference  model  used  10  meter  grid  spacing  and  a  time  step  of  1 
millisecond.  These  grid  and  time  step  dimensions  are  well  within  the  accepted  guidelines 
for  numerical  stability  and  dispersion  requirements  of  the  grid  formulation  (Kelly  et  al., 
1976;  Virieux,  1986).  The  length  of  the  run  (6  seconds  or  6000  timesteps)  was  long 
enough  for  steady  state  to  be  reached. 

The  quantity  of  interest  from  the  finite  difference  runs  is  the  scattered  pressure  field. 
Since  the  total  field  is  the  default  output  of  the  method,  the  unperturbed  source  field  must 
be  subtracted  to  leave  the  scattered  field.  This  was  done  by  computing  an  acoustic 
halfspace  model  (no  cylindrical  scatterer)  with  receivers  in  the  same  positions  as  in  the 
scattering  model.  Then,  the  total  and  unperturbed  finite  difference  fields  are  both  known 
and  the  scattered  field  can  be  obtained  by  a  simple  subtraction.  The  time  series  of  the 
scattered  field  for  the  CW  model  are  shown  in  figure  4.3.  This  is  the  scattered  pressure  in 
the  water  at  a  distance  of  1.4  kilometers  from  the  center  of  the  cylinder.  The  time  series  are 
plotted  against  azimuths  from  180  degrees  (backscattering  direction)  to  360  degrees  (full 
-orward  scattering  direction).  Scattering  angles  from  1  to  179  degrees  have  symmetric 
counterparts  in  the  series  presented. 

The  time  series  demonstrate  a  number  of  features  unique  to  the  scattering  from  a 
CW  source.  First,  the  transition  to  steady  state  takes  place  about  2.0  seconds  after  the 
beginning  of  the  scattered  signal  (or  about  4.25  seconds  into  the  time  series).  Most  of  the 
ime  series  go  through  some  amplitude  fluctuations  until  about  4.25  seconds  when  the 
implitude  levels  out  Since  the  CW  source  is  monochromatic,  the  time  series  also  contain  a 
single  frequency  for  all  azimuths  and  times.  Subsequent  to  the  onset  of  steady  state,  the 
only  quantity  in  the  time  series  which  is  variable  is  the  pressure  amplitude  as  a  function  of 
azimuth.  Because  of  the  elastic  nature  of  the  cylinder  and  the  existence  of  circumferential 
waves  (in  the  full  wave  solution  of  the  finite  difference  method)  there  is  no  geometric 
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Figure  4.3.  Time  series  of  scattered  pressure  from  the  10  Hz  CW  source  for  circular 
receivers  between  180  and  360  degrees  (see  figure  4.2)  at  a  distance  of  1.4  km.  from  the 
cylinder  center.  These  data  show  a  strong  dependence  of  scattering  strength  on  azimuth. 
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Figure  4.4.  Analytical  (solid  line)  vs.  finite  difference  (dashed  line)  results  of  scattered 
pressure  strength  vs.  azimuth  for  the  10  Hz  CW  source  model.  Cylinder  radius  is  0.160 
km. 
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shadow  zone  in  the  forward  direction.  In  fact,  the  strongest  signal  occurs  in  forward 
directions  and  smaller  amplitudes  occur  at  a  number  of  azimuths  within  the  figure.  This  is 
also  predicted  by  the  normal  mode  analytical  solution. 

Both  the  analytical  and  finite  difference  results  for  scattering  from  an  elastic 
infinite  cylinder  with  the  10  Hz  CW  source  are  shown  in  figure  4.4.  The  analytical 
solution  was  calculated  by  solving  equation  15  for  azimuths  between  1  and  360  degrees. 
The  value  for  the  finite  difference  solution  are  the  maximum  pressure  at  each  azimuth  taken 
after  steady  state  was  reached.  Again,  because  of  symmetry,  only  half  of  the  problem  was 
solved  with  finite  differences  and  plotted  in  figure  4.4.  The  striking  feature  of  the  pattern 
for  CW  scattering  from  a  cylinder  is  the  fact  that  energy  is  scattered  primarily  into  distinct 
azimuths.  Cylinders  which  are  very  large  or  very  small  with  respect  to  the  seismic 
wavelength  do  not  show  this  type  of  pattern  as  dramatically.  The  lobes  in  the  scattering 
pattern  are  really  due  to  an  interference  pattern  created  by  the  different  wave  types  of  the 
system.  Body  wave  such  as  the  direct  plane  wave  source,  waves  transmitted  through  the 
elastic  cylinder  (as  well  as  whispering  gallery  modes),  and  reflected  waves  interfere  with 
each  other  and  also  with  circumferential  waves  of  both  Franz  and  Rayleigh  types  to  create 
the  patterns  seen  in  the  figure.  Franz  waves  are  due  solely  to  Huygens  principle  and  exist 
for  perfectly  rigid  cylinders  but  do  not  exist  for  a  flat  or  gently  sloping  interface  (without 
shadow  zones).  Rayleigh-type  circumferential  waves  exist  only  with  elastic  media  and 
degenerate  to  true  Rayleigh  waves  for  flat  interfaces. 

The  agreement  between  the  analytical  and  finite  difference  solutions  is  quite  good 
considering  the  coarseness  of  the  grid  used  and  the  fact  that  a  rectangular  rather  than 
cylindrical  coordinate  system  is  used  for  the  finite  difference  formulation.  A  curved 
surface,  such  as  the  edge  of  the  cylinder,  defined  on  a  rectangular  grid  must  necessarily  be 
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stepwise.  Previous  work  (chapter  3)  has  shown  that  scattering  does  occur  from  these  small 
steps  which  will  affect  the  interference  pattern  seen  in  figure  4.4. 

The  finite  difference  method  also  does  some  inherent  averaging  over  space  (and 
time  )  in  calculating  spatial  derivatives.  This  averaging  can  blur  the  actual  location  of 
boundaries.  The  interference  pattern  seen  in  a  plot  of  scattered  pressure  vs.  azimuth  is  very 
sensitive  to  the  diameter  of  the  cylinder  for  a  given  wavelength  of  energy.  Analytical 
scattering  patterns  for  cylinders  with  radii  of  one  grid  step  more  and  less  (0.17  and  0.15 
km.)  and  one  half  grid  step  more  and  less  (0. 165  and  0.155  km)  that  that  used  for  the  finite 
difference  model  are  shown  in  figure  4.5.  Small  changes  in  diameter  (or  ka)  have  a 
dramatic  influence  on  the  scattering  pattern  and  the  averaging  done  by  the  method  probably 
does  account  for  at  least  some  of  the  discrepancies  between  the  analytical  and  finite 
difference  solutions,  especially  for  azimuths  near  the  backscattered  direction 
(90<azimuth<27 0) . 

It  is  extremely  difficult,  if  not  impossible,  to  pick  contributions  from  the  individual 
wave  types  in  the  time  series  or  scattering  patterns.  A  few  of  the  individual  contributions 
can,  however,  be  seen  in  the  snapshots  of  compressional  energy  on  the  grid.  Figure  4.6  is 
the  snapshot  of  compressional  and  shear  scattered  (unperturbed  source  field  removed) 
energy  on  the  grid  4.0  seconds  after  introduction  of  the  CW  source  onto  the  top  of  the  grid. 
Normalized  compressional  and  shear  energy  was  calculated  by  using  the  divergence  and 
curl  of  displacements  output  by  the  finite  difference  formulation  (see  Dougherty  and 
Stephen,  1988  (chapter  2)).  The  scattered  field  was  obtained  from  the  total  field  by  the 
same  procedure  used  for  the  time  series  of  subtracting  the  homogeneous  halfspace  solution 
from  the  total  field  to  leave  the  scattered  fieid. 
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Figure  4.5.  Analytical  (solid  line)  vs.  finite  difference  (dashed  line)  results  of  scattered 
pressure  strength  vs.  azimuth  for  the  10  Hz  CW  source  model.  Cylinder  radii  of  0.15  km 
(4.5a.),  0.155  km.  (4.5b.),  0.165  km.  (4.5c.),  and  0.170  km.  (4.5d.)  for  the  analytical 
patterns.  The  finite  difference  results  are  for  a  cylinder  radius  of  0. 16  km.  The  analytical 
pattern  is  especially  sensitive  to  cylinder  radius  in  the  backscattered  direction. 
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Figure  4.6.  Snapshot  of  compressional  (top)  and  shear  (bottom)  scattered  energy  for  the 
10  Hz  CW  source  finite  difference  model,  ^-is  snapshot  was  taken  4.0  seconds  after  the 
source  was  turned  on.  The  cylinder  is  located  on  the  axis  of  symmetry  on  the  left  hand 
boundary  of  the  griu.  Compressional  energy  in  the  water  is  scattered  only  into  a  few 
discrete  directions.  Shear  energy  only  appears  inside  of  the  elastic  cylinder. 
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The  main  scattering  features  which  can  be  seen  in  the  snapshots  of  figure  4.6  are 
the  interference  pattern  and  the  energy  inside  the  cylinder.  The  interference  panem  in  the 
fluid  manifests  itself  as  energy  scattered  into  certain  discrete  azimuths  with  the  main  lobe 
being  in  the  forward  direction  (no  'shadow  zone').  Also,  at  distances  greater  than  about  3 
cylinder  diameters,  the  scattered  field  does  not  change  significantly  with  increased  distance 
from  the  cylinder  center  (along  a  particular  azimuth).  Two  clear  contributors  to  the 
scattered  interference  pattern  in  the  fluid  can  be  seen  within  the  elastic  cylinder  and  along  its 
boundary.  Body  waves  transmitted  into  the  cylinder  can  be  seen  inside  the  elastic  solid 
(refer  to  figure  4.2  for  cylinder  boundaries).  Although  much  of  the  transmitted  energy 
probably  is  reverberated  as  'whispering  gallery'  modes,  direct  and  vertically  travelling 
compressional  and  shear  waves  can  be  seen  in  the  upper  and  lower  plots  respectively. 
These  waves  interfere  directly  with  circumferential  waves  both  inside  and  outside  of  the 
cylinder  boundary.  Elastic  components  of  the  interface  waves  can  be  seen  along  the  inside 
of  the  cylinder  boundary  especially  in  the  shear  energy  plot  (lower  part  of  figure  4.6). 
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PULSE  SOURCE 


Another  method  which  can  be  used  to  examine  scattering  from  an  elastic  cylinder  is 
to  use  a  broad  band  pulse  source  rather  than  a  monochromatic  CW  source.  Although  the 
scattering  patterns  presented  in  figures  4.4  and  4.5  are  useful,  they  only  provide  the 
response  for  a  single  frequency.  Another  quantity,  Foo,  the  scattering  form  function, 

provides  the  strength  of  scattering  in  any  direction  for  all  frequencies  present  in  the  broad 
band  source.  Each  monochromatic  CW  solution  contains  one  value  of  Foo  at  a  particular 
value  of  ka.  By  using  a  broad  band  pulse  source,  a  range  of  values  of  Foo  can  be 
established  in  a  single  experiment  or  finite  difference  run.  In  this  way,  a  pattern  of  the 
influence  of  frequency  (or  conversely,  diameter  of  the  cylinder  for  a  single  frequency)  on 
the  scattering  strength  can  be  determined.  Since  it  is  known  that  scattering  strength, 
especially  in  the  backscattering  direction,  is  very  sensitive  to  the  diameter  of  the  cylinder 
(see  figures  4.4  and  4.5)  perhaps  a  better  test  of  the  finite  difference  method  is  to  examine 
the  pattern  of  the  backscattering  form  function  rather  than  try  to  precisely  match  a  particular 
monochromatic  response. 

A  time  domain  solution  to  the  transient  source  problem  can  also  give  added  insight 
to  the  mechanisms  of  scattering  which  result  in  the  interference  patterns  of  the 
monochromatic  source  problem.  The  short  duration  of  the  pulse  allows  the  wavefield  to 
separate  into  discrete  components  after  initial  interaction  with  the  cylinder.  Time  series  for 
both  the  analytical  and  finite  difference  solutions  can  be  obtained  and  compared  for  the 
existence  and  location  of  arrivals  in  time.  Successive  snapshots  of  wavefronts  show  quite 
clearly  how  the  incident  energy  is  partitioned  between  direct,  transmitted,  and 
circumferential  phases.  Another  advantage  of  using  the  pulse  source  is  that  steady  state 
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need  not  be  obtained  for  the  analysis  (although  sufficient  time  must  be  run  to  include  all  or 
the  significant  arrivals). 

Theoretically,  the  use  of  a  broadband  pulse  source  to  determine  the  amplitude  of  the 
backscattering  form  function  is  straight  forward  and  follows  the  work  of  (Dardy  et  al., 
1977;  Neubauer,  1986).  The  plane  wave  pressure  source  used  for  the  finite  difference 
models  is  in  the  shape  of  the  second  derivative  of  a  Gaussian  curve  and  is  given  by  (Kelly 
et  al.,  1976;  Stephen,  1985); 

(16)  pj  (t)  =  A  [  4  ^  (  3  t  -  2  ^  t3 )  exp  (  f  )  ] 

and  the  Fourier  transform  is  given  by; 

(17)  gj  (ka)  =  j  pi  (t)  exp  ( -jwt  ) 

The  parameter  Z,  in  equation  (16)  determines  the  peak  frequency  and  bandwidth  of  the 
source  wavelet.  Both  the  source  time  function  and  its  Fourier  transform  are  given  in  figure 
4.7.  This  wavelet  has  a  peak  frequency  of  10  Hz  and  upper  and  lower  halfpower 
frequencies  of  13.5  and  6.8  Hz  respectively.  The  peak,  upper,  and  lower  frequencies 
correspond  to  ka  values  of  6.7, 9.0,  and  4.5  respec,ively. 

The  form  function  for  a  given  receiver  azimuth  relates  incident  to  scattered  pressure 
for  a  range  of  ka  values.  From  equation  (15),  the  scattered  pressure  in  the  far  field  for  a 
given  azimuth  is; 
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and  the  form  function  for  the  backscattered  direction  is  given  by; 
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The  backscattered  pressure  has  the  Fourier  transform; 
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so  that  the  amplitude  of  the  backscattered  form  function  can  be  expressed  as; 
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The  finite  difference  method  uses  a  pressure  source  function  given  by  equation  (10)  and 
calculates  scattered  pressure  (equation  1 8).  After  calculation  of  the  Fourier  transforms  of 
these  two  pressures  (equations  17  and  20)  the  amplitude  of  the  backscattered  form  function 
with  respect  to  ka  can  be  determined  by  equation  21.  The  analytical  solution  to  scattering 
of  the  source  function  given  in  equation  i 6  by  an  aluminum  cylinder  can  be  obtained  by 
equation  (15).  Calculation  of  the  backscattered  form  function  for  this  case  by  using 
equation  (21)  agrees  exacdy  with  the  analytical  curve  for  an  aluminum  cylinder  seen  in 
Dardy  et  al  (1977,  figure  3). 
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Figure  4.7.  Frequency  content  (4.7a.)  and  pressure  time  signal  (4.7b.)  of  the  Gaussian 
pulse  source  used  for  the  transient  problem. 
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Figure  4.8.  Grid  geometry  used  for  the  finite  difference  solution.  Absorbing  boundaries 
are  used  on  the  top  and  right  sides  of  the  grid.  The  left  hand  side  is  an  axis  of  symmetry. 
Energy  travelling  downward  is  attenuated  in  the  stipled  region  by  using  the  telegraph 
equation.  The  cylinder  is  placed  in  the  center  of  the  grid  for  the  transient  problem  to 
facilitate  detection  of  circumferential  waves  travelling  around  the  cylinder  boundaries. 
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The  cylinder  model  was  rerun  with  a  few  slight  modifications  for  the  pulse  source. 
The  geometry  for  these  models  is  shown  in  figure  4.8.  Although  more  computer  memory 
and  time  is  needed,  the  cylinder  was  placed  at  the  center  of  the  grid  in  order  to  observe  the 
predicted  circumferential  waves  travelling  around  the  water-solid  interface.  Receivers  are 
again  placed  every  degree  in  circles  centered  about  the  cylinder  center.  The  plane  wave, 
broadband  pulse  source  (figure  4.7)  is  introduced  along  the  top  of  the  grid  and  has  only 
vertical  particle  motion. 

Because  of  the  apparent  problem  of  scattering  from  the  coarse  boundary  observed 
in  the  CW  model,  pulse  source  models  with  both  5  (PULSE5)  and  10  (PULSE  10)  meter 
grid  spacing  were  computed.  The  finer  5  meter  grid  spacing  causes  the  cylinder  edge  to  be 
smoother  when  discretized  onto  the  finite  difference  grid.  This  allows  for  a  reduction  in  the 
edge  diffraction  problem  of  the  coarser  models.  The  overall  dimensions  of  the  10  and  5 
meter  models  are  300x400  grid  steps  x6000  times  steps  and  600x800  grid  steps  x  12000 
time  steps,  respectively. 

Snapshots  of  pressure  determined  by  the  finite  difference  method  are  shown  in 
figure  4.9  for  model  PULSE5.  The  time  domain  solution  provided  by  the  method  gives 
added  insight  into  the  wavefronts  which  contribute  to  the  problem  which  is  not  available 
with  the  CW  source  models.  The  partitioning  of  energy  from  the  pulse  is  readily  visible  in 
the  snapshots.  Snapshots  in  figure  4.9  are  shown  for  every  0.15  seconds  starting  0.5 
seconds  after  the  plane  wave  source  is  introduced  along  the  top  of  the  grid.  In  figure  4.9a, 
the  plane  wave  source  has  just  reached  the  cylinder.  Elastic  wave  velocities  within  the 
cylinder  (p-wave  velocity  of  2.6  km/sec,  s-  wave  velocity  of  1 .5  km/sec)  cause  the  portion 
of  the  wavefront  transmitted  into  the  cylinder  to  move  ahead  of  the  source  wave  and 
increase  in  wavelength. 
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The  scattered  wavefield  becomes  much  more  complicated  after  the  source  wave  has 
passed  completely  through  the  region  of  the  cylinder  (figures  4.9b-f).  The  circular 
wavefronts  which  emanate  from  the  cylinder  in  later  time  steps  are  complex  combinations 
of  reflected,  transmitted,  and  circumferential  waves  predicted  by  various  ray  and  acoustic 
theoretical  studies. 

In  the  backscattered  direction  (towards  the  top  of  each  snapshot  in  figure  4.9),  the 
strongest  and  most  obvious  scattered  arrival  is  the  reflected  wave  seen  in  all  but  the  first 
snapshot  of  figure  4.9.  Energy  which  is  transmitted  once  through  the  cylinder  makes  up 
the  strongest  arrivals  in  the  forward  scattering  directions.  Most  of  the  important 
wavefronts  which  combine  in  the  CW  case  to  form  the  complex  interference  patterns  of 
figures  4.4  and  4.5  are  present  in  figure  4.9d.  The  reflected  (R)  and  transmitted  (T)  waves 
are  results  of  simple  ray  theoretical  interactions  with  the  cylind'^r  boundaries.  Fast  (FI)  and 
slow  (SI)  interface  or  circumferential  waves  are  due  to  the  curved  nature  of  the  interface. 
The  slow  interface  wave  is  similar  to  the  Franz  wave  present  in  the  rigid  cylinder  case 
(Doolittle  and  Uberall,  1966;  Frisk  et  al.,  1975;  Karal  and  Keller,  1964).  In  this  elastic 
ca^e,  however,  it  does  have  an  elastic  wave  component  inside  the  cylinder.  This 
circumferential  wave  has  a  cardioid  shape  in  the  fluid  portion  of  the  model  and  shows  up 
best  in  figures  4.9d  and  4.9e.  Whispering  gallery  modes  and  the  interference  of  fast 
interface  waves  make  up  the  high  energy  pulse  travelling  in  the  backscattered  direction  just 
inside  of  the  reflected  wavefront  (seen  in  figure  4.9c,d,f)-  Other  energy  in  the  snapshots  is 
due  to  continuous  leaking  of  whispering  gallery  modes  back  into  the  fluid  and  scattering  of 
these  modes  and  circumferential  waves  from  the  comers  present  on  the  discretized  cylinder 
boundary. 
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Figure  4.9.  Snapshots  of  compressional  energy  in  the  finite  difference  grid  of  the  pulse 
source  model  PULSE5.  See  figure  4.8  for  the  location  of  the  cylinder  in  the  center  of  the 
grid.  Cylinder  radius  is  0.16  km.  with  grid  spacing  of  5  meters.  The  snapshot  time 
increment  is  0.15  seconds.  In  the  first  frame  (4.9a.),  the  source  has  just  illuminated  the 
cylinder  causing  a  backscattered  reflection  and  a  transmitted  wavefront  which  is  slightly 
ahead  of  the  source  plane  wave  because  of  the  faster  P-wave  velocity  in  the  cylinder.  0.75 
seconds  later  (4.9f.),  the  scattered  field  is  comprised  of  a  complex  combination  of  reflected 
and  diffracted  waves,  whispering  gallery  reverberations,  and  circumferential  waves 
(cardioid- shaped  wavefront  in  figure  4.9f.). 
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All  of  these  arrivals  can  also  be  seen  in  the  time  series  for  circular  receivers  located 
1.4  km.from  the  center  of  the  cylinder  (figure  4. 10).  Time  series  traces  can  also  be 
obtained  analytically  for  receivers  at  the  same  discrete  points  and  are  shown  in  figure  4. 10. 
Theoretically,  snapshots  such  as  those  produced  by  the  finite  difference  method  could  also 
be  calculated  analytically  but  would  take  a  very  large  amount  of  computer  time  given 
present  computational  resources.  Analytical  time  series  are  created  for  each  angle  by 
multiplying  the  scattering  form  function  (equation  19  at  any  azimuth)  and  source  power 
(equation  17)  in  the  frequency  domain  for  all  frequencies  present  in  the  pulse  source.  An 
inverse  Fourier  transform  of  equation  (20)  (for  a  given  azimuth)  then  gives  the  scattered 
pressure  time  series  at  the  receiver, 
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Equation  (22)  was  solved  for  integer  angles  up  to  360  degrees  and  plotted  in  figure 
4. 10a.  Time  series  arrivals  are  marked  with  the  same  notation  as  in  the  snapshots  of  figure 
4.9.  Circumferential  waves  cross  in  the  backscattered  direction  (azimuth=  1 80  degrees)  as 
they  travel  around  the  cylinder  boundary.  The  results  for  PULSE5  and  PULSE10  are  also 
shown  in  figure  4.10.  Two  obvious  numerical  features  are  present  in  the  finite  difference 
results  and  can  be  ignored.  These  are  the  high  amplitude  source  arrivals  in  the  lower  rught 
and  left  hand  sides  of  figures  4.10b  and  4.10c  (the  analytical  solution  only  contains  the 
scattered  pressure,  not  the  total  field)  and  the  axis  of  symmetry  arrivals  in  the  upper  left 
hand  side  of  figures  4. 10b  and  4. 10c.  All  of  the  arrivals  expected  from  the  analytical  time 
series  are  also  present  in  the  finite  difference  results.  PULSE5  has  good  agreement  but 
there  are  discontinuities  in  the  arrivals  and  some  additional  incoherent  arrivals  due  to  edge 
scattering. 
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As  mentioned  earlier,  the  backscattered  form  function  is  often  used  to  describe 
scattering  from  a  cylinder.  This  is  obtained  analytically  by  using  equation  21.  The  finite 
difference  value  for  this  is  derived  from  the  signal  of  the  receiver  at  the  backscattered 
direction  with  the  source  pulse  removed.  The  procedure  used  is  essentially  the  inverse  of 
the  procedure  used  for  creating  the  analytical  time  series  but  with  the  finite  difference 
scattered  pressure  in  the  left  hand  side  of  equation  22.  These  functions  for  the  analytical 
and  PULSE5  models  are  shown  in  figure  4.1 1  for  ka  values  within  the  band  of  the  source 
frequencies.  Shown  in  the  figure  are  analytical  (solid  line)  traces  assuming  cylinder  radii  of 
0.16  km,  0.165  km.  and  0.17  km.  (figures  4.11a,  4.11b,  4.11c,  respectively)  compared 
with  the  finite  difference  solution  for  a  cylinder  radius  of  0.16  km.  (dashed  lines  in  figures 
4. 1  la,  4. 1  lb,  4. 1  lc).  The  shapes  of  the  curves  for  the  finite  difference  model  do  not 
change  noticeably  but  the  different  cylinder  radii  values  have  the  effect  of  a  phase  shift  on 
the  curves.  The  best  match  for  the  curves  is  found  for  a  cylinder  radius  of  0. 17  km. 
Apparently,  the  finite  difference  source  wave  is  'seeing'  a  slightly  larger  cylinder  than  is 
actually  being  used.  In  all  cases,  the  absolute  level  of  the  backscattered  form  function 
curve  is  comparable  between  the  analytical  and  finite  difference  solutions.  Considering  the 
sensitivity  of  the  backscattered  pressure  to  cylinder  radius  seen  in  figure  4.4  and  4.5,  the 
finite  difference  approximation  to  the  cylinder  does  an  excellent  job  reproducing  scattering 
levels  and  effects. 
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Figure  4.10.  Analytical  (4.10a.)  and  finite  difference  (PULSE5  in  4.10b.  and  PULSE10 
in  4. 10c.)  pressure  time  series  vs.  azimuth  for  the  transient  source  problem.  The  finite 
difference  results  contain  two  numerical  features  which  do  not  appear  in  the  analytical 
solution.  These  are  the  return  from  the  axis  of  symmetry  (curved  arrival  in  the  upper  left 
hand  portion  of  figure  4.10b.  and  4.10c.)  and  the  direct  source  arrival  (strong  arrivals  in 
the  lower  right  and  left  hand  portions  of  4. 10b.  and  4. 10c.). 
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Figure  4. 1 1 .  Backscatter  form  function  for  analytical  results  (solid  line'*  and  model 
PULSE5  (dashed  line).  Cylinder  radii  used  to  plot  PULSE5  results  (see  equation  21)  is 
0.160  km.,  b.  0.165  km.,  and  c.  0.170  km. 
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BACKSCATTER  FORM  FUNCTION 


CONCLUSIONS 


1.  The  problem  of  acoustic  plane  wave  scattering  from  an  elastic  cylinder  is  a  stringent  test 
of  the  finite  difference  method  in  rectangular  coordinates.  The  formulation  used  here, 
along  with  the  grid  parameters,  is  the  same  as  that  used  for  earth  models  and  provides  good 
agreement  to  analytical  results  of  scattering  from  cylinders. 

2.  The  method  is  computationally  intensive  but  the  solution  obtained  is  a  complete  one 
which  qualitatively  matches  that  given  by  the  normal  mode  solution. 

3.  The  offset  grid  used  in  the  finite  difference  solution  provides  the  analytically  correct 
symmetric  solution  for  the  plane  wave  sources. 

4.  These  results  verify  the  observation  by  other  authors  that  the  problem  of  backscattering 
from  a  cylinder  is  very  sensitive  to  cylinder  diameter  (or  conversely,  to  frequency).  The 
finite  difference  solution  appears  to  be  seeing  a  slightly  larger  cylinder  that  that  actually 
defined  on  the  grid. 

5.  Rough  comers  caused  by  the  definition  of  the  curved  cylinder  boundary  on  the  square 
finite  difference  grid  complicate  the  solution.  Decreasing  the  size  of  the  grid  step  lessened 
this  problem  for  the  pulse  source  models  but  finer  grids  are  probably  necessary  for  an 
'exact'  match  to  the  analytical  solution.  This  problem  is  heightened  by  the  magnitude  of  the 
velocity  contrast  at  the  cylinder  boundary  and  is  probably  not  as  important  for  the  earth 
models  with  much  smaller  velocity  heterogeneity  contrasts  and  gradational  (rather  that 
sharp)  scatterer  boundaries. 
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Conclusions 


The  subject  of  this  work  is  the  scattering  of  seismo/acoustic  energy  from  lateral 
heterogeneities  in  the  upper  oceanic  crust.  For  simplicity,  most  seismic  studies  in  the  past 
have  assumed  lateral  homogeneity  even  though  obvious  scattering  effects  are  seen  in  almost 
all  seismic  field  data.  Another  common  theme  throughout  the  chapters  of  this  thesis  is  that 
the  scatterers  studied  all  have  spatial  dimensions  on  the  same  order  of  magnitude  as  the 
seismo/acoustic  wavelength.  In  this  realm  of  scattering  problems,  scattering  reaches  a 
maximum  and  full  wave,  laterally  heterogeneous  modeling  schemes  are  essential.  The 
finite  difference  method  of  calculating  wave  propagation  through  heterogeneous  elastic 
media  was  used  for  the  various  studies  comprising  this  work.  Results  from  the  modeling 
demonstrate  the  importance  of  considering  lateral  heterogeneity  as  a  scattering  mechanism. 
Also,  the  work  of  chapters  three  and  four  point  out  some  restrictions  of  the  method  for  use 
with  the  rough  seafloor  scattering  problem.  A  brief  summary  of  the  results  from  the 
individual  chapters  is  presented  here. 

Chapter  one  dealt  with  the  modeling  of  a  diffracted  arrival  seen  in  a 
particular  line  of  ocean  bottom  hydrophone  data  from  the  Rivera  Ocean  Seismic 
Experiment.  This  arrival  is  referred  to  as  a  'refraction  branch  diffraction'  because  of  its 
location  in  time  and  space  after  the  refracted  arrival  and  to  distinguish  it  from  a  near  normal 
incidence  diffraction  hyperbola.  Finite  difference  models  of  deterministic  seafloor  features 
demonstrated  that  this  arrival  could  be  due  to  diffraction  from  any  of  a  number  of 
topographic  features  in  the  area.  Energy  scattered  from  the  seafloor  scatterer  travels 
through  the  crust  back  to  the  ocean  bottom  receiver.  Other  arrivals,  such  as  converted 
shear  waves,  interference  waves  and  the  pseudo-Rayleigh  wave  occur  in  the  models  but  not 
in  the  original  field  data.  This  is  probably  due  to  the  use  of  an  inappropriate  Poisson's  ratio 
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for  the  area  (no  shear  conversion  for  higher  Poisson's  ratio)  and  the  presence  of  more 
complicated  topography  in  the  field  area  (which  would  attenuate  the  interference  waves). 

The  problem  of  incoherent  scattering  from  continuous  randomly  heterogeneous 
upper  oceanic  crust  was  investigated  in  chapter  two.  Although  these  models  dealt  with  flat 
seafloors  and  the  strength  of  the  heterogeneity  may  be  debatable,  this  type  of  crust  is 
probably  much  more  realistic  that  a  strictly  laterally  homogeneous  one.  For  random 
velocity  perturbations  with  a  Gaussian  distribution,  there  is  a  general  trend  of  increasing 
random  scatter  as  ka  (k=wavenumber,  a=correlation  length)  approaches  one.  Secondary 
deterministic  scattering  from  the  heterogeneities  increases  as  ka  increases  above  one  (larger 
heterogeneities).  Self-similar  distributions  cause  scattering  with  characteristics  of  a  range 
of  ka  values. 

Heterogeneities  in  the  models  presented  in  chapter  two  act  as  secondary  sources  for 
Stoneley  waves  along  the  water-solid  interface.  This  mechanism  may  be  one  explanation 
for  the  generation  of  seafloor  noise  which  may  also  propagate  partly  in  the  form  of  interface 
waves.  The  final  major  result  from  chapter  two  is  that  propagation  through  continuously 
heterogeneous  crust  affects  the  correlation  of  primary  arrival  waveforms.  In  general, 
decorrelation  reaches  a  maximum  as  ka  approaches  one. 

The  effects  of  scattering  from  seafloor  topography  can  be  very  similar  to  those  from 
lateral  heterogeneities.  The  purpose  of  the  work  in  chapter  3  was  to  investigate  scattering 
from  more  complicated  topographic  profiles  than  the  isolated  deterministic  scatterers  of 
chapter  one.  A  more  important  result  from  chapter  3  was  the  evaluation  of  the  use  of  the 
finite  difference  method  with  a  rectangular  grid  for  rough  seafloor  problems.  The  problems 
encountered  in  using  this  method  involve  the  discretization  of  a  curved  interface  on  a 
rectangular  grid.  There  is  significant  scattering  of  energy  from  the  stair-step 
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microroughness  superimposed  on  the  larger  scale  topography.  Although  this 
microroughness  should  not  be  present  to  solve  the  sinusoidal  seafloor  problem,  the 
relatively  large  amount  of  energy  scattered  from  these  small  features  (1/15  wavelength) 
does  point  out  the  need  to  consider  small  scale  topography  in  the  scattering  problem. 
Because  some  of  this  energy  also  travels  as  interface  waves,  this  could  be  another 
mechanism  for  the  generation  of  seafloor  noise. 

Despite  the  problem  of  scattering  from  the  microroughness,  many  'real'  effects  of 
propagation  of  the  sinusoidal  topography  were  also  observed  in  the  models.  Travel  time 
anomalies,  multiple  compressional  and  shear  diving  waves,  as  well  as  some  strong 
backscattered  arrivals  are  due  to  the  larger  scale  topography.  Steep  topography  allows 
energy,  especially  shear  energy,  to  enter  the  seafloor  even  at  great  ranges.  Ray  theoretical 
shadow  zones  are  not  present  because  of  Franz-type  waves  diffracted  into  areas  where  the 
grazing  angle  is  less  than  zero. 

Chapter  4  involved  the  use  of  the  time  domain  finite  difference  method  to  solve  the 
problem  of  acoustic  wave  scattering  from  an  infinite  elastic  cylinder.  Scattering 
mechanisms  and  grid  considerations  seen  in  this  chapter  are  closely  tied  to  those  seen  in  the 
results  of  chapter  three.  The  same  finite  difference  formulation  and  grid  used  in  chapters 
two  and  three  was  used  unmodified  for  the  cylinder  problem  in  chapter  4.  This  formulation 
provided  good  agreement  to  the  analytical  solution  of  the  problem.  The  results  obtained 
here  agree  with  the  observations  of  other  authors  that  this  problem  is  very  sensitive  to 
cylinder  diameter.  The  finite  difference  solution  appears  to  be  seeing  a  slightly  larger 
cylinder  than  that  actually  defined  on  the  grid. 

Discretization  of  the  round  cylinder  cross-section  on  the  rectangular  finite  difference 
grid  caused  the  same  type  of  problems  as  seen  in  chapter  three  with  sloping  seafloors. 
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These  problems  were  lessened  somewhat  in  the  models  of  both  chapters  by  decreasing  the 
grid  step  size  but  an  even  finer  grid  would  be  necessary  to  more  accurately  represent 
smoothly  varying  interfaces.  The  magnitude  of  the  velocity  contrast  between  the  solid  and 
the  surrounding  acoustic  medium  heightens  this  problem  of  scattering  from  rough  comers 
(caused  by  the  discretization).  This  is  probably  not  an  issue  for  the  earth  models  in  chapter 
two  with  dramatically  smaller  velocity  contrasts  and  gradational  scatterer  boundaries. 

An  important  point  brought  out  by  the  results  of  chapters  three  and  four  is  that  the 
complexity  of  the  model  must  be  matched  to  the  appropriate  modeling  technique.  The 
'simple'  models  presented  in  these  chapters  are  probably  better  solved  using  deformed 
grids,  different  coordinate  systems,  or  analytical  methods.  However,  the  objective  of  this 
work  was  to  evaluate  the  use  of  the  finite  difference  method  with  a  rectangular  grid  for 
rough  seafloor  problems.  It  was  found  that  extremely  fine  grids  must  be  used  in  order  to 
accurately  represent  the  smoothly  varying  curved  interfaces  of  the  sinusoidal  seafloor  and 
cylinder  models.  The  power  of  the  finite  difference  method  lies  in  its  ability  to  handle 
models  with  media  too  complex  to  be  solved  analytically.  The  models  in  this  study  are  a 
starting  point  for  more  complex  models  and  help  to  point  out  some  of  the  restrictions  of  the 
method.  If  the  specification  of  the  medium  must  be  at  a  very  fine  scale  for  detailed  earth 
structure,  then  a  very  fine  finite  difference  grid  is  necessary  and  appropriate.  The  finite 
difference  method  using  a  rectangular  grid  will  be  an  important  tool  for  determining  just 
how  detailed  the  seafloor  must  be  sampled  deterministically  to  accurately  represent 
interaction  with  a  seismo/acoustic  wavefield.  It  may  turn  out  that  interaction  with  smaller 
features  (with  respect  to  energy  wavelength)  can  be  represented  as  a  stochastic  component 
superimposed  on  the  deterministic  effects  of  larger  scale  topography. 
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Common  themes  occur  throughout  the  chapters  of  this  work.  First,  since  the  sizes 
of  the  heterogeneities  in  all  of  the  models  have  at  least  one  dimension  which  is  on  the  order 
of  magnitude  of  the  seismic  wavelength,  scattering  effects  are  strong  and  concern  the  wave¬ 
like  characteristics  of  the  energy.  Scattering  effects  from  lateral  volume  heterogeneities  and 
small  topographic  features  (such  as  the  grid- imposed  steps  of  chapters  three  and  four)  are 
similar  and  both  types  of  heterogeneities  can  act  as  secondary  sources  for  Stoneley  waves 
along  the  water-solid  interface.  Deterministic  effects  such  as  travel  time  and  amplitude 
anomalies  are  stronger  for  the  rough  seafloor  models  because  of  the  sharp  impedance 
contrast  and  size  of  the  large  scale  topography.  These  effects  would  be  less  for  sediment 
covered  or  'softer'  seafloors  and  smaller  topographic  features.  Franz  waves  and  the  lack  of 
shadow  zones  due  to  the  interaction  of  acoustic  waves  with  the  curved  interface  of  a 
cylinder  can  also  be  seen  in  the  sinusoidal  seafloor  models. 

The  work  presented  here  ties  in  well  with  current  areas  of  interest  in  marine 
seismology.  There  is  much  interest  presently  in  the  generation  and  propagation  of  seafloor 
noise.  It  has  been  shown  here  that  scattering  from  topographic  and  volume  heterogeneities 
can  be  an  important  mechanism  for  coupling  of  body  waves  (both  in  the  water  and  in  the 
upper  crust)  into  interface  waves  along  the  seafloor.  Areas  of  future  work  in  this  area 
include  determining  the  relative  importance  of  volume  vs.  surface  scattering  mechanisms 
and  the  influence  of  bottom  elastic  parameters  on  the  generation  of  secondary  Stoneley 
waves.  Closely  spaced  array  experiments  should  be  able  to  identify  modes  of  propagation 
of  the  noise. 

Also  of  current  interest  is  the  fine  scale  characterization  of  the  upper  oceanic  crust. 
The  models  shown  here  lead  to  the  hypothesis  that  there  is  information  regarding  the  fine 
scale  structure  of  the  oceanic  crust  to  be  found  in  the  secondary  features  of  seismograms. 
Inversion  of  seismic  data  for  heterogeneity  length  scales  should  be  possible  after  further 
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evaluation  of  model  results  such  as  those  presented  here.  This  may,  however,  involve  the 
collection  of  seismic  data  with  very  closely  spaced  receivers  in  order  to  properly  identify 
phases.  The  ultimate  goal  of  obtaining  knowledge  of  fine  scale  crustal  structure  is  to 
understand  mechanisms  of  crustal  formation  at  the  mid-ocean  ridges  and  the  processes  of 
evolution  of  the  crust  with  age. 

Reverberation  of  energy  from  the  seafloor  is  another  area  of  research  activity  where 
these  results  are  relevant.  Not  only  have  the  models  in  chapter  three  shown  the  importance 
of  small  scale  (with  respect  to  wavelength)  topographic  features  in  the  scattering  problem, 
they  have  also  demonstrated  an  important  mechanism  for  allowing  energy  to  enter  the  upper 
crust  at  large  ranges  (beyond  the  flat  seafloor  critical  range).  It  is  currently  not  understood 
whether  'reverberated'  energy  is  due  to  scattering  from  seafloor  topography  or  volume 
heterogeneities.  The  type  of  modeling  done  here  will  help  to  determine  the  importance  of 
considering  elastic  parameters  of  the  bottom  as  well  as  seafloor  topography  in  the  acoustic 
reverberation  problem.  Computer  resources  are  becoming  available  which  will  allow  very 
fine  grids  to  be  used  in  rough  seafloor  models.  These  grids  will  be  used  to  determine  the 
influence  of  fine  scale  topography  on  the  scattered  field,  that  is,  whether  or  not  the  small 
scale  interactions  must  be  modelled  deterministically  or  whether  they  can  be  included  as  a 
stochastic  component  superimposed  on  larger  scale  deterministic  effects. 
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